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We study the quantum transport through a strongly interacting Anderson quantum dot at finite 
bias voltage and magnetic field at zero-temperature using the real-time renormalization group (RT- 
RG) in the framework of a kinetic (generalized master) equation for the reduced density operator. 
To this end we further develop the general finite temperature real-time transport formalism by 
introducing field superoperators that obey fermionic statistics. This direct second quantization 
in Liouville-Fock space strongly simplifies the construction of operators and superoperators which 
transform irreducibly under the Anderson-model symmetry transformations. The fermionic field 
superoperators naturally arise from the univalence (fermion-parity) superselection rule for the total 
system of quantum dot plus reservoirs. Expressed in these field superoperators, the causal structure 
of the perturbation theory for the effective time-evolution superoperator-kernel becomes explicit. 
Using the constraints of the causal structure, we construct a parametrization of the exact effective 
time-evolution kernel for which we analytically find the eigenvectors and eigenvalues in terms of a 
minimal set of only 30 independent coefficients. The causal structure also implies the existence of a 
fermion-parity protected eigenvector of the exact Liouvillian, explaining a recently reported result 
on adiabatic driving [Phys. Rev. B 85, 075301 (2012)] and generalizing it to arbitrary order in 
the tunnel coupling T. Furthermore, in the wide-band limit the causal representation exponentially 
reduces the number of diagrams for the time-evolution kernel. The remaining diagrams can be 
identified simply by their topology and are manifestly independent of the energy cut-off term- 
by-term. By an exact reformulation of this series we integrate out all infinite-temperature effects, 
obtaining an expansion targeting only the non-trivial, finite-temperature corrections, and the exactly 
conserved transport current follows directly from the time-evolution kernel. From this series the 
previously formulated, RT-RG equations are obtained naturally. We perform a complete 1- plus 
2- loop RG analysis at finite voltage and magnetic field, while systematically accounting for the 
dependence of all renormalized quantities on both the quantum dot and reservoir frequencies. Using 
the second quantization in Liouville-space and symmetry restrictions we obtain analytical RT-RG 
equations with an efficient numerical solution and we extensively study the model parameter space, 
excluding the Kondo regime where the 1 plus 2-loop approach is obviously invalid. The incorporated 
renormalization effects result in an enhancement of the inelastic cotunneling peak, even at a voltage 
~ magnetic field ~ tunnel coupling Y. Moreover, we find a tunnel-induced non-linearity of the 
stability diagrams (Coulomb diamonds) at finite voltage, both in the SET and ICT regime. 

PACS numbers: 73.63.Kv, 05.10.Cc, 03.65.Yz, 05.60.Gg 



I. INTRODUCTION 

Non-linear transport spectroscopy of nanoscalc sys- 
tems is a key technique in modern day physics. Although 
well-understood in linear transport (equilibrium)^, see 3 - 
for a review, the theoretical description of the non- 
equilibrium transport regime, especially at low temper- 
ature, remains challenging and the recent progress in 
this direction has led to an improved understanding of 
quantum transport through strongly interacting systems. 
For a recent comparative review see^. Recent years 
have seen the appearance of many fully numerical ap- 
proaches, such as the scattering-state time-dependent 
numerical renormalization group (TD-NRG)^ relying 
on the discretization parameter approximation 6 -, time- 
dependent density matrix RG (TD-DMRG) 7 - - — , iterative 
path integrals (IPSI)^, numerically-exact influence func- 
tional path integrals (INFPI)ii^, quantum Monte Carlo 
(QMC) in combination with Nakajima-Zwanzig projec- 



tion technique^ or with imaginary-time formalism^ and 
diagrammatic Monte Carlo (diagMC)^. Partially ana- 
lytical approaches involve the non-crossing approxima- 
tion (NCA) 16 i 17 ! equations of motion for Green func- 
tions^, Bethe-Ansat a 19 ' 20 and the flow-equation ap- 
proach 2 ^. Fully analytical approaches include Keldysh 
perturbation theory to high orders in [ / 22 i 23 or dual- 
fermion superperturbation theory 2 ^. Finally, perturba- 
tive renormalization group studies have mostly started 
from the Kondo-model mapping of the Anderson model, 
cither working with Kcldysh-Grcen functions^ or the re- 
duced density operator approac h 26 ! 27 . 

The density operator approach has a long history in 
various fields in physics and chemistry. In the context of 
nanoscale transport, it is a natural starting point for the 
description of systems with large interaction energies in 
the high temperature - weak coupling limit, U 3> T ^>Y. 
It can be systematically extended to include high-order 
tunneling processes using the real-time diagrammatic^ 
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or Nakajima-Zwanzig approach^ 9 - technique, in partic- 
ular when combined with Liouville space techniques 3 ^. 
However, at low temperature this approach becomes 
problematic due to effect of high-energy contributions 
which rcnormalize the effective low-energy physics. In 
this paper we show how this approach can be extended 
to this regime. In particular, we show that much of 
the insightful structure of generalized master / kinetic 
equation approach at high temperature is preserved in 
this regime. We formulate a real-time renormalization 
group approach that naturally connects to the general- 
ized quantum master or kinetic equation approach. 

For this purpose we study the simplest possible, bench- 
mark model of an interacting quantum dot coupled to 
metal electrodes, the Anderson model. We explicitly 
set up a general approach to deal with the non-linear 
transport at very low temperature making use of real- 
time perturbation theory and a renormalization group 
formulated based on this. We proceed analytically as far 
as possible, making this a technically challenging task. 
In fact, the problem is not really manageable without a 
new physical approach for dealing with superoperators. 
The development of this Liouville- Fock space approach 
is a central topic of this paper. Our approach differs 
by both construction and the scope of the application 
from the previous formulations^—. The most close to 
the our approach is one introduced by Prosen^i, which 
was used to calculate steady states of quadratic effec- 
tive Louvillians. Here we extend it to the reservoirs with 
continuous fields as well and further develop it to sim- 
plify the microscopic derivation of effective Liouvillians 
for non-quadratic problems. The approach introduced 
by Schmuta 3 -^ and used by other author s 33 ' 34 differs from 
our approach in principal details which are essential for 
our application. 

A large part of the paper is devoted to this new ap- 
proach since it allows for general physical insights into the 
problem and is crucial for overcoming of technical diffi- 
culties in setting up RT-RG. We will illustrate its poten- 
tial in application to RT-RG calculation of the non-linear 
transport at zero temperature. As result this paper by 
necessity extensive. Further motivation for this is that 
several general physical insights into the the real-time 
approach haven not been pointed out, although this ap- 
proach has been developed for some time and have found 
widespread use, see2£ for a review. To include these it is 
required that the approach is set up from scratch, paying 
special attention to: (i) causal structure and the related 
Kcldysh rotation (ii) Liouville space formulation, in par- 
ticular a 2nd quantization for superoperators (iii) spin 
and charge rotation symmetry (iv) the infinite tempera- 
ture limit as a reference point for both the 2nd quantiza- 
tion technique that we develop, as well as the perturba- 
tion theory. Only by fully exploiting these the simplest 
application of RT-RG to the Anderson model becomes 
practically possible. Clearly, these developments are best 
presented coherently in the context of the application to 
the RT-RG for which it is absolutely crucial. 



To indicate the impact of these developments for the 
Anderson model study, we note that the simplest approx- 
imation which includes the exact result for the U = 
limit requires an infinite series of diagrams in the stan- 
dard RT-PT. Using the RT-RG this result is recovered 
only when performing a 1- and 2-loop analysis for the 
effective Liouvillian and including 1-loop vertex correc- 
tions. We emphasize that when applied to the interact- 
ing case, this incorporates renormalization effects from 
strong tunneling, while neglecting spin-fluctuation pro- 
cesses relevant only in the Kondo regime, which enter 
only in a 3-loop RG analysis (the latter has been ad- 
dressed previously based on a Kondo- model mapping^ 7 -) . 
Naively formulating these equations leads to hundreds of 
non-linear coupled integro-differential RG equations for 
frequency-dependent coupling functions. The central re- 
sult of this paper is the derivation of 30 coupled differen- 
tial coupling functions systematically incorporating the 
leading frequency dependence which includes the U = 
limit exactly. On the way we derived several exact results 
of general importance. This makes an efficient numerical 
implementation possible and allows experimentally rele- 
vant stability diagrams to be calculated from wide ranges 
of parameters in the non-linear zero-temperature regime 
(excluding the narrow Kondo regime). 

Many of the results can be extended to generic models 
involving local interactions (multi-orbital Anderson type 
mod els) with bilinear tunnel coupling to reservoirs. Sec 
Ref. l35| for a recent study of non-local interactions using 
the RT-RG. 

The paper is organized in two main parts as follows. 
In the first part, starting in Sec. pi)) , we formulate the 
model and directly change to a Liouville space descrip- 
tion and develop the kinetic equation approach for the 
stationary QD density operator. We formulate the per- 
turbation series for the effective Liouvillian L(z) appear- 
ing in this equation using a new causal representation of 
field superoperators G with fermionic statistics. We em- 
phasize that this formulation of the perturbation theory 
- although equivalent to previous formulation s 29 ' 30 ' 36 - 
leads to many simplifications beyond the application of 
interest here, and therefore warrants a proper, extensive 
discussion. Several of these results have already found 
application 3 ^, and even provide insights into and gener- 
alization of recent interesting predictions 3 ^. Moreover, a 
rcnormalized perturbation theory, which takes the infi- 
nite temperature limit as a formal reference point, sug- 
gests itself. It also connects in a natural way to the renor- 
malization group approach while preserving much of its 
general perturbative structure. We show how the calcu- 
lation of the current requires little additional calculation 
and prove that in our non-linear approach the linear cur- 
rent vanishes at zero bias, which is not obvious from the 
general structure of the theory. 

In the second main part of the paper, starting in 
Sec. (|IIip. the explicit 1- and 2-loop RT-RG equations are 
derived accounting for the energy dependence of both the 
Liouvillian and vertices due to the finite non-equilibrium 
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transport voltage. The non-interacting current is shown 
to arise naturally as an exact result already in 1-loop 
RT-RG, implying that all 2-loop corrections to this ob- 
servable arise from the strong local Coulomb interaction. 
We find, however, that for U = non-zero 2-loop terms 
exist which are relevant when one is interested in, e.g., 
the density matrix (and not just the current). The non- 
trivial frequency dependence of 2-loop equations is sys- 
tematically accounted for in powers of the renormalizcd 
dimensionless coupling supcropcrators G, resulting in an 
effective RG equation for the effective Liouvillian only, 
which accounts for vertex renormalization corrections. 

This simplification enables the detailed numerical 
study in Sec. (|IVj) of the zero-temperature non-linear 
transport in all regimes excluding the Kondo regime 
of low applied voltage and magnetic field. The impor- 
tance of accounting for both 1- and 2-loop corrections 
as well as the Matsubara axes is demonstrated numer- 
ically. Finally, in Sec. (|IVj) we show that the tunnel- 
induced renormalization effects incorporated in our 1 plus 
2 loop approach enhance the inelastic cotunneling reso- 
nance at finite magnetic field and voltage and generate 
non-linearities of the SET stability diagrams (Coulomb 
diamonds). 



II. MODEL AND REAL-TIME TRANSPORT 
THEORY 



In the continuum limit the reservoirs are described by 
the density of states v r {uS) — J^k H UJ ~ e r.k+t l r) an d we go 
to the energy representation of the fermionic operators, 

a<j.r(u) = 1 V" a ay r. k S(uj - e r ,k + Mr), (5) 

with the anti-commutation relations: 

[a a ^(Lo),a\, r ,(uj')] + = 5 a ^>5 T y5(u - u') (6) 
[a a , r {w),a a < y(u)')} + = 0. (7) 

where we denote (anti)commutators by [A, B]± = AB ± 
BA. Thus we have for the reservoir Hamiltonian: 

H R = Yj [ dLJ ^ + Mr)4,r(">K,rM (8) 
<7,r 

In contrast to e r< k the energy oj is the electron energy 
relative to /i r , i.e., the reference energy depends which 
reservoir is considered. 

The junctions connecting the dot and reservoirs are 
modeled by the tunneling Hamiltonian 

V = Y / V r (9) 

r 

^ = J2 I <W^rH (i r M4, r (wK + h.C.) (10) 



A. Anderson model 

In this section we introduce the model and our compact 
notation which is crucial to the Liouvillc space formula- 
tion of the theory. The simplest model Hamiltonian of a 
QD which takes into account Coulomb interaction effects 
involves just a single orbital 



H = en + BS Z + Un^n^ 



(1) 



Here e denotes the energy of the orbital with occupa- 
tion operators n = n a , n a = d\d a . The index a = ± 
corresponds to spin f, I and S z = \ o~n a is the z com- 
ponent of the spin vector operator S = ^2 aa i \(T a ,a'd\d a i 
along the external magnetic field B = Be z (in units 
where g[iB = 1) and <r is the vector of Pauli matrices. 
The dot is attached to electrodes which are treated as 
free electron reservoirs: 



H 1 



E 



(2) 



The Hamiltonian of the total system is denoted by 



H tot = H + H R + V 



(11) 



We assume t r (uj) to be real and introduce the spectral 
density 



T r (oj) — 27T ■v r (uj)\t r {uj)\ 2 
and rescaled field operators: 



/ r r(w) 



(12) 



(13) 



To make the notation more compact we introduce an 
additional particle-hole index^ 



& t,rH' 7 l z 

b a , r (u), n -- 
~ 1 d a , n = - 



(14) 
(15) 



where the reservoir index r = ± corresponds to L, R, k is 
the orbital index and a the spin index, quantized along 
the z-axis. The reservoir electron number and spin can 
be decomposed into n R = n r and = ^ r s r where 



E a &,r,k a >?,r,k 
a.k 



,r,k 



(3) 
(4) 



Throughout the paper we will denote the inverse value 
of a two- valued index with a bar. A prominent example, 
is given by 



V = -v 



(16) 



We combine all indices into a multiindex variable written 
as a number: 



l = r),o-,r,u; l = *n,a,r,w. 



;i7) 
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where by way of exception the bar denotes inversion of 
the particle-hole index only. Then b\ = b„ ar (ui) and 
b\ = &_^ |CT)r (a;) and the various independent anticommu- 
tation relations are compactly summarized by 

[d u d 2 }+ = 5 1 2 (18) 

[61,62]+ = ^*ia (19) 

Z7T 

where Ti = T r (uj). The interaction then simply reads 

V = hdi (20) 

where we implicitly assume summation over all discrete 
parts of the multi-index 1 (i.e. 77, c, r) and integration 
over its continuous part {uj). We usually omit subscripts 
in the multiindcx components 77, a, r, uj if we have the only 
multiindcx. If we have more than one multiindex we dis- 
tinguish their components by corresponding subscripts: 
1 = ?7i,cri,ri,cJi, 2 = 7]2,<J2,r 2 ,uj2- Importantly, opera- 
tors (and, below also superoperators) of the dot and the 
reservoirs can treated as if they commute (rather than 
anticommute)££. 

The reservoirs are assumed to be at thermal equilib- 
rium with temperature T, each described by its own 
grand-canonical density operator, 

P R = X[p r , P r = ^e-^ Hr -^ nT \ (21) 

r 

where z r = Tre~T( Hr ~^" r ). We assume that a symmet- 

r 

ric bias is applied to the electrodes, i.e., pl,r = ±V/2. 
We note the key property 

b x p R = e^/VV (22) 

In Eq. (fl"2j) the density of states varies on the energy 
scale of the bandwidth D which we assume to be much 
larger than any other energy scale in the problem. In this 
wide-band limit we can assume T r (uj) to be energy inde- 
pendent and cut off all reservoir energy integrals (uj) at 
the scale D. The detailed energy dependence of T r (uj) at 
high energies is not crucially important for the results^. 
In the actual applications in Sec. (|IV[) we will assume for 
simplicity that the tunnel couplings are symmetric, i.e. 

= Tr = r and consider the low temperature limit i.e. 
T < U, V, T r by setting T = 0. The results of the present 
section and much of Sec. . however, do not depend 
on these assumptions unless explicitly indicated. 

For a non-zero magnetic field B and finite Coulomb 
interaction U the total system possesses two locally and 
globally conserved quantities that will play an important 
role. Locally on the dot and the reservoirs the charge 
and spin component along the magnetic field on the dot 
are conserved, 

[H,n}- = [H,S z ]-=0 (23) 
[H R , n R ]- = [H R , S R ]- = (24) 



These conservation laws extend to the total charge, 
N tot = n + n R , and spin, 5* ot = S z + S R , since the 
interaction V commutes with these operators: 

[H tot ,N tot ]- = (25) 

[jjtot )jS .tot]_ = Q ( 2 g) 

B. Density operator and diagram rules 

The purpose of this section is twofold. Firstly, we 
briefly review the real-time approach to the calculation 
of the stationary reduced density operator, introducing 
the central quantities S(z), the self-energy supcropcra- 
tor, and L(z), the effective Liouvillian, and their per- 
turbativc expansions in vertex superoperators G. Sec- 
ondly, we introduce a "causal" representation of the per- 
turbation theory which allows for a new compact formu- 
lation and derivation of the diagrammatic rules for the 
self-energy E. Moreover, many general physical insights 
become explicitly clear in this representation. In partic- 
ular, starting from this new formulation, the possibility 
of a two-stage real-time renormalization group (RT-RG) 
that will be set up in Sec. piB4|) and Sec. (|IIIj) arises 
quite naturally. 

1. Stationary density operator 

In order to find the QD stationary state we need to con- 
sider the evolution of the total system density-operator. 
It evolves according to the Liouville- von Neumann equa- 
tion: 

d t p tot {t) = -i [H tot ,p tot (t)] _ = -iL tot p tot (t) (27) 

with the superoperator Liouvillian L tot — [H tot , •] _ . Su- 
peroperators are linear transformations of operators and 
throughout the paper we let • indicate the operator on 
which a superoperator acts (if needed). Explicit matrix 
representations of superoperators are only required for 
the QD part and will be discussed later on in Sec. (jll C[) 
The initial state of the total system at the initial time 
to is assumed to be the direct product of the dot density 
matrix and the equilibrium density matrices (|2ip of the 
electrodes: 

p tot (t )=p(t )p R (28) 

We will discuss some properties of p(to) further below. 
The formal solution of Eq. ([27]) is: 

p tot (t) = e -iH^(t-t ) p tot( to yH*°\t-t ) (2g) 

= e- iLtot ^p tot (to) (30) 

and the reduced density matrix of the dot is obtained by 
integrating out of reservoirs degrees of freedom: 

p{t) = Trp(t) = Tr (e"^ ^ p(t ) p R ) (31) 
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We now decompose L tot = L+L R +L V , with L = [#,•]_, 
L R = [H R , •]_ and set up the perturbation series in the 
tunnel coupling L v = [V,»]_ ~ vT 7 . It is then more 
convenient^ to use the Laplace transform of the dot re- 
duced density matrix for Imz > 0: 



when we define field superoperators as follows: 



p(z) = J dte iz ^ to) pit) 

to 

= Tr / 



(32) 



P tot (to) (33) 



R \ z-L-L R -L v > 



We will refer to z as the dot frequency. We expand the 
resolvent in L v , resulting in a geometric series with terms 
of the form 



Tr 

R 



1 



L — L R 



I 



L R 



P to \t . 



(34) 



The average over the reservoirs can now be calculated 
directly using a Wick theorem for field superoperators 
(see Eq. ([501) an d App.rS}. Collecting irreducible contrac- 
tions into the self-energy supcropcrator £(z) (see below, 
Sec. pi B 3J) ) the perturbation series can be resummed: 



(35) 



z- L(z)' 

where we have introduced the effective Liouville operator 

L(z) = L + Y,(z) (36) 

To keep the notation to a minimum we distinguish this 
quantity from the "bare" dot Liouvillian L = [H, •] 
by simply appending the dependence on the frequency 
z. L(z) completely determines the time-evolution of 
the reduced density operator. A key idea exploited 
both in the perturbation theory and the renormaliza- 
tion group is that one is free to modify L and as 
long as their sum remains equal to L(z). The equa- 
tion determining the stationary density matrix p = 

limt_t -Kx> pif) = lim2^.;o(~ iz)pi z ) is now obtained by 
multiplying Eq. ((35)) by — iziz — Liz)) and taking z — > iO: 



L(i0)p= 



(37) 



Before deriving the perturbation series for S(z) in 
Sec. piB 3[) . we first introduce a convenient represen- 
tation of the field superoperators. 



2. Causal representation of fermionic field superoperators 

Below we integrate out explicitly the reservoir degrees 
of freedom while keeping track of those of the QD. To 
facilitate this the tunnel coupling supcropcrator L v = 
[V, •] _ should be written as a convenient product of the 
dot and the reservoir superoperators: inserting Eq. (|20[) 
we have 



V 



P 



(38) 



di», 



(39) 
(40) 



The superscript p = ± keeps track of whether a field 
operator acts from the left or right and is referred to 
as Keldysh index by analogy to the Green-function and 
functional integral techniques. In Eq. (f3"5)) we implicitly 
sum over p, in addition to the multiindex 1. A crucial 
difference to the formulation of Ref. [26| is that we intro- 
duced a harmless additional supcropcrator 



P 



into Eq. (|38|). where 



P 



N tot 


= [iV tot ,.]_ 




L n 


= [n,»]-, 






= \n R ,.U 





(41) 

(42) 
(43) 

(44) 

arc the superoperators associated with the total, QD 
and reservoir electron numbers respectively. Clearly for 
p = +1 this factor is trivially equal to 1. However, it 
may seem at first sight that for p = — 1 this is not the 
case: when applied to a projector of states of the total 
system this superoperator counts the relative parity of 

the fermion numbers N and N : (-1) L |iVA)(iV'A'| = 
±\NX)(N'X'\ for N — N' = even / odd, where A, A' de- 
notes further quantum numbers. However, in all calcu- 
lations we can assume that the total-system state opera- 
tors on which it acts have even parity, since odd fermion- 
parity components of states can neither be measured by 
any physical operator nor be prepared using physical evo- 
lutions. This is referred to as the fermion-parity supers- 
election rul o 39 i 40 . 

Since the fermion parity plays an important role in 
what follows but is often not mentioned or used explic- 
itly in density operator approaches it warrants some dis- 
cussion, in particular since odd-fcrmion parity operators 
do appear in the renormalization group approach. Phys- 
ical Hamiltonians and observablcs (and their correspond- 
ing superoperators) always contain only products of even 
numbers of fermionic operators. This implies that only 
the part of the density operator p tot it) with even fermion 
parity can enter into the calculation of any physical ob- 
servable (A) it) = Tr Ap tot it). This even part of p tot (t) 
is generated solely from the even parity part of /9 tot (to) 
at earlier times since by the same token the parity of 
the total fermion number is conserved during time evo- 
lution. Therefore only the even fermion-parity part of 
the initial state p tot (io) can contribute to an observable 
and one may set any odd-fermion part of any density ma- 
trix equal to zero. As a result, we can take in Eq. (|38[) 

JV tot 

p • = 1», even for p = — 1. For the factorized ini- 
tial state Eq. (|28p that we assumed here, this implies that 
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p(to) must be assumed to be of even fermion parity, since operators, 
p R also has even fermion parity (L n n p F 



[n R ,p R ] 



0, 

which follows from Eq. (j21|l). 

The useful implications of this fermion-parity conser- 
vation become clear when performing a linear transfor- 
mation of the field supcroperators with respect to their 
Keldysh indices p. The naively chosen field superopera- 
tors ([39| and (|40|) have the disadvantage that they com- 
mute or anticommute depending on the Keldysh index 
p: 



(45) 
(46) 



This complicates many calculations as noted, e.g., 
in Ref.lU However, the factorization of the total fermion 
parity into a dot and reservoir factor in (|41[) naturally 
suggests a transformation of the field operators. By ab- 
sorbing the fermion parity superoperators of each sub- 
system into new field superoperators, 



P LnR J 



they obey definite anticommutation relations 



[QIQIY 



-+ =p8 pp ,5 1 2 



(47) 
(48) 

(49) 
(50) 



This allows one to prove the Wick theorem directly for 
the operators J using simple algebra^, avoiding the need 
to carefully keep track of sign factors as done in Ref. I2H . 
However, the non- vanishing anticommutators still de- 
pend on the Keldysh index p on the right hand side. 
Moreover, the fields have no simple Hermitian supercon- 
jugation relation (see App. |B]). This can be avoided by a 
rotation of the QD fields 



G q 



q = 

p 



(51) 

and a contravariant rotation (cf. Sec. (jll B 3[l ) of the 
reservoir fields 



T 9 



■h 



f{ q 



. (52) 

To facilitate later discussions, we introduce both an index 
q = ± as well as ~ and ~ symbols to distinguish the new 
field components. Now the anticommutation relations 
are completely analogous to those of the usual fcrmionic 



r<i ni 
u l7 u^, 



J + 

Gi,G v ] =0 
J + 

^ [G 1 ,Gv} + =0 



(53) 



whereas in the reservoirs we incorporate the coupling into 
the normalization factor: 



T q J q 



£1 
2tt 



8q,q>Sl,V 



Jl, Jv 
Jl, J\i 

[Ji,Jv 



+ 2tt 
= 



(54) 

This second transformation is known as the Keldysh rota- 
tion^ and is, e.g., applied to fermionic fields represented 
by Grassman-numbers in functional integral theories^ 
or in Green-function formalism^ 2 -. Here we find that it 
also considerably simplifies the real-time transport the- 
ory in which not all degrees of freedom can be integrated 
out, in contrast to the cited approaches. We note that 
a transformation similar to Eqs. (|4"Tl) - (f4"5)) which also re- 
sults in the usual anti-commutation and conjugation re- 
lations for the fermionic superoperators superoperators 
was introduced in Ref. 13214341 however, without perform- 
ing of the Kcldysh-rotation (see also discussion in ApplDl 
). This transformation, however, is less convenient in gen- 
eral since it does not reveal a general structure of the 
fermionic superoperators which is important for the our 
applications and which we will relate to causality below. 

We therefore will refer to (|5"Tj) - ([52"j) as the causal rep- 
resentation in Keldysh space and the index q = ± as the 
causal index. Again we denote its inverse by q = —q. 

For the dot field superoperators (G q ) this representa- 
tion was already introduced in Ref. l26l but it was not ex- 
ploited in the context of the perturbation theory, where 
it provides many useful additional simplifications that we 
now address. First of all, the property 



Tr(G.) 

D 







(55) 



ensures the probability conservation on the dot during 
the time-evolution, see Eq. (|69[) . It is extremely impor- 
tant for the formulation of the RT-RG in Sec. and is 
preserved during renormalization. We note that proba- 
bility conservation in the Liouville approach corresponds 
to the normalization conservation for the partition func- 
tion in the path-integral approach^. Eq. (|5"5j) follows 
directly from the definition of the causal representation: 
0, and, using the cyclic property of the 



using TrL™ 

D 

trace, Tr&? 

D 

fore 



Tr(G.) 



Trdi» is independent of p 



El r P 



±. Thcre- 



(56) 
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Secondly, the fields G\ and Gj are related by Hcrmitian 
conjugation of supcropcrators: 

Gi = {Gi?. (57) 

For a superoperator S the Hermitian conjugate is de- 
fined by Tr (A* SB) = Tr((5U)tS) where are A and B 

are arbitrary operators, see also Sec. pi C 1[) and App. [B] 
This indicates that the conjugate fields Gi and Gj arc 
similar to usual creation and annihilation operators, re- 
spectively. In Sec. pi C|) we exploit this "2nd quantiza- 
tion" in Liouville space to construct a convenient basis 
the QD Liouville space that includes the left and right 
zero super-eigenvectors of both vertices. This consider- 
ably simplifies the matrix representations of many supcr- 
opcrators required in perturbation theory and the RG. 
An immediate consequence of the probability conserva- 
tion property ([57jl is that G has the unit operator 1 as a 
right zero super eigenvector: 

Gil = (58) 

The above properties of these causal field supcropcr- 
ators acting on the QD Liouville space provide explicit 
insight into important physical issues that are otherwise 
not obvious in the general form of the perturbation the- 
ory, e.g., the wide-band limit, the infinite temperature 
limit, and the energy cut-off D dependence which we 
discuss in Sec. piB 3p . Introducing the corresponding 
causal representation for the reservoir field supcropcra- 
tors also yields several simplifications which we now dis- 
cuss. 

a. Wick theorem First, we note that due to the local 
interactions on the QD Eq. ([53"]) cannot be used to for- 
mulate a Wick theorem for the vertices G. In contrast to 
this, for the reservoir field superoperators this is possible 
due to the relation 

J lP R = t£mh(r ll Lu 1 /2T)J 1 p R (59) 

This result follows from Eq. (|22p by writing it in superop- 
erator notation, p R = e VlLJl ^ T J*^ p R , and then ap- 
plying the transformations (gSJl and flSjj). With Eq. (|59"|) 
we can algebraically prove the Wick theorem for the su- 
peroperators Jf in close analogy to the usual case of 
fermionic operators (see App. |A"]) : the average of a prod- 
uct equals the product of pair contractions summed over 
all contractions of pairs (ik), 

Tr(jf ...Jtp R ) =EIl(- 1 ) P ( J f 4")R (60) 

i<k (ik) 

with the usual fermionic sign (— l) p of the permutation 
P that disentangles the contractions. 

b. Causal structure The number of possible pair- 
contractions appearing on the right hand side of Eq. 

is strongly reduced. Applying TrJi» to Eq. (|59")l written 

for the operator J 2 we obtain a relation between two of 
the possible four pair contractions: 

(JiJ2)r = tanh(r ?2 w 2 /2T)(J 1 J 2 ) iJ (61) 



This is actually a statement of the equilibrium 
fluctuation-dissipation theorem for each reservoir sepa- 
rately, see, e.g., Ref. 0- In close analogy to Eq. (|55|) one 
proves 

Tr(Ji«) = (62) 

This implies in particular, (JiJ%)r = which is related 
to the well-known fact that of the four reservoir Green's 
functions only three are independent^. This is analogous 
to the absence of so called "quantum-quantum" contrac- 
tions (or "classical-classical" term in the Kcldysh action) 
in the Kcldysh functional integral approach and has been 
referred to as the "causal" structure of Green functions 
in that approach^. However, Eq. (|()2"j) entails an addi- 
tional simplification: (J\ J 2 )_r = 0. This is particular to 
our real-time superoperator approach, and is most ex- 
plicitly related to causality. To see this, note that we 
keep track of the left and right action of operators us- 
ing superoperator notation and that these superopera- 
tors inherit their ordering in the Laplace representation 
from their forward time-ordering in the time-evolution, 
cf. Eq. (f3"4"|) . We never have to introduce a fictitious 
backward time-propagation as in the Kcldysh technique. 
Therefore no advanced reservoir Green functions can ap- 
pear in our theory, which is expressed by (JiJ 2 )h = 0. 

c. Energy and temperature dependence The only 
reservoir correlation functions that can appear in the 
Wick expansion (|60|) for the causal fields have a sim- 
pler energy (lo) and temperature (T) dependence than in 
the representations (|3"9" ]) -(|3D ]l and (|17 ]) -([4"8" ]) . This corre- 
lates with the physical information that these functions 
incorporate, the retarded function^ 

r 2 

71,2(772^2) := (J\J2)r = 7^—^1,2, (63) 

Z7T 

and Kcldysh function 

- - r 2 

7i )2 (?7 2 w 2 ) := (JiJ 2 )r = — tanh(r7 2 w 2 /2T)5 1 2 . (64) 

Z7T 

In the causal representation the superoperator ordering 
explicitly shows that the retarded contraction 7 con- 
tains no information about the distribution function of 
the reservoirs and is therefore temperature independent. 
Eq. (|63|) follows directly from the anticommutation rela- 
tion ([51| and the property (jol?)) . and is indeed indepen- 
dent of the reservoir density operator p R . In contrast, 
for the Keldysh contraction, one first needs to use the 
property (|59[) specific to the equilibrium state of the non- 
interacting reservoirs, before Eq. (|5"4")) and Eq. (|6"2"|) can 
be applied. The above shows that the representation (f52|) 
reflects most explicitly the causal structure of the pertur- 
bation theory, motivating its denotation. It thereby au- 
tomatically achieves the decomposition of the reservoir 
Fermi distribution function into its symmetric (trivial) 
and antisymmetric (non-trivial) part with respect to the 
energy u> that was introduced in Ref. [26l . where it is cru- 
cial in setting up the RT-RG. Below we show that in the 
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causal representation for the QD fields the perturbation 
theory also drastically simplifies and that in this new for- 
mulation the RT-RG appears quite naturally. 



3. Perturbation series and diagram rules 

a. Diagram rules When defining the causal repre- 
sentation (Eq. (jSTj) and Eq. (|52")l ) we performed opposite 
Keldysh rotations for the dot and the reservoirs. This 
is motivated by the form of the tunneling Liouvillian 
Eq. (|38p , which in the causal representation can be writ- 
ten compactly as 



L v = OlJi + GJi = J2 G l J l 



(65) 



q=± 



We can now integrate out the reservoir degrees of freedom 
in each term of the expansion (|34[) in the way discussed 
in Ref. [26| by commuting all reservoir operator to the 
right side using the relation 



J q L 



resulting for the order m contributions to S(^) in 



1 



z-L 
1 



Gf 



1 



-G 



<n 



1 



X n 



z-Xt-L A z-X 2 -L 
-G q " 



(66) 



(67) 



p{t )Jr{jfJf..J!tp R ) 

Z — Li H 



Here X t = £ k<i Xk where the summation runs over the 
reservoir frequencies x k — Vk(^k + Mfc) 01 the G^ fc (</?*) 
originally standing to the left to the resolvent i. Applying 
the Wick theorem (|60|) we can represent the terms dia- 
grammatically by propagators connecting vertices that 
are contracted in pairs by lines with frequencies Xk- The 
irreducible parts of these diagrams, i.e., those parts which 
cannot be cut without hitting at least one reservoir con- 
traction, arc collected into the irreducible kernel or self- 
energy superoperator S(z). In Fig. [T] we show the dia- 
grams for S(z) to 1- and 2-loop order. Leaving the sum 
over all possible configurations of pair contractions, all 
indices and all orders m implicit we can write: 



i 



(68) 



z-Xi-L 



-G q2 C? 9 ™- 1 - 



z — X m _i — L 



-G q " 



Here 7 denotes that the function 7 (7) , given by Eq. ([64jl 
(Eq. (|63|1 ). should be written for pair contraction con- 
necting a G vertex on the left with a G (G) vertex on the 
right. By Eq. (|65p it is the earliest vertex (rightmost) 
that decides the type of contraction, i.e., its indices ap- 
pear as the argument of the contraction function. Here 
Xi is now the sum over the frequencies of all reservoir 
contractions which go over i-th resolvent (since contrac- 
tions that start and end to the left and to the right cancel 
out). 



Importantly, since on the left we always have G (i.e. 
G qi with qi = — ), the property (|55j) of the causal field 
of type G is seen to guarantee the conservation of prob- 
ability: 



TrS(z) = 

D 



(69) 



Another general property restricts the frequency 
dependence^ 



KY,(z)K = -£(-; 



(70) 



Here K = K^ 1 is the anti-linear superoperator that 
effects the Hermitian conjugation of an operator, see 
App. [G] This guarantees through Eq. (J3SJ) - dSHJ) that 
the reduced density operator remains Hermitian during 
the time-evolution, i.e., p(z) = p{—z*). This property 
is guaranteed by the conjugation properties (|120[) and 
KLK^ 1 = —L. The above causal representation of 
the diagrammatic perturbation theory is very useful in 
general and is extended and applied extensively to time- 
dependent problems in Ref. 1371 . We now discuss the main 
advantages, which will be important for setting up the 
RG and the construction of a convenient supervector ba- 
sis in Sec. (|HC|) . 

b. Wide-band limit First of all, the number of con- 
tributing terms is strongly reduced since the trivial (7) 
and non-trivial (7) energy dependence of the contractions 
is automatically separated which in other representations 
would have to be done separately^: Terms which do not 
contribute in the large bandwidth limit can be identified 
diagrammatically, using their diagram topology as illus- 
trated in Fig. [TJ all diagrams where one or more vertices 
are enclosed between contracted vertices Gi and Gj give 
contributions of order T/D <C 1 or smaller and can be 
neglected. This can be proved by careful examination 
of the poles appearing when closing all integrals in the 
complex upper half plane, see Ref. 1371 or argued in the 
time-representationj21 As a result 7 contractions can only 
occur inside 7 contractions in diagrams with more than 1 
loop, see Fig. [1] This results in an exponential reduction 
in the number of contributing terms. 

This feature also naturally suggests a starting point for 
a two stage RG approach. First, all terms which contain 
retarded contractions 7 can be integrated out explicitly 
by a one-step diagram resummation, leading to a renor- 
malization of L discussed in Sec. piB4[) . The remaining 
diagrams will then contain only G dot supcropcrators 
with non-trivial contractions 7 which require a second, 
continuous RG. This approach will be worked out in de- 
tail in Sec. (|HT|) . 

c. Infinite-temperature limit The choice of the su- 
pervector basis is simplified very much by noting that 
the superoperator structure of T,(z) is strongly restricted 
in the causal representation: in the wide-band limit (see 
above) there is only one diagram that starts with G on 
the right and ends with G on the left. We denote this 
special 1-loop diagram by S. Importantly, all other dia- 
grams start and end with a G vertex. The uj integral for 
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Gj Gg G2 Gi Gj Gi G2 Gi 



+ 



Gi G5 G2 G 



1 K?2 Gj G2 Gi 




C?2 ^1 ^2 (jl 



T/D < 1 



- T 2 /D 2 1 



FIG. 1: 1- and 2-loop contributions to the self-energy ker- 
nel E. The full black line denotes the free dot propagation, 
the arrow indicates the ordering of superoperators ( "late goes 
left") in Laplace space that is "inherited" from the time- 
ordering. The curves denote the contraction of dot super- 
operators. Here the line connects the dot superoperators 
whose corresponding reservoir superoperator are contracted: 

(GJ... GJ) R = (JJ) R G ...G — > G...G . This is in con- 
trast to standard techniques where a contraction line connects 
contracted operators themselves. The full blue line denotes a 
Keldysh contraction 7 of two G type vertices (cf. Eq. (|64[) ~) 
whereas the dashed blue contraction corresponds to the re- 
tarded contraction 7 of a G with an earlier G vertex to the 
right (cf. Eq. (|63[0 . Diagrams with vanishing contractions, 
i.e., where G is contracted to the right (cf. Eq. (|62[0 . are 
not drawn. Furthermore, the 2 loop diagrams in the light 
and dark gray boxes contain a 7 contraction enclosing k = 1 
and respectively k — 2 other vertices and therefore scale with 
F k /D k and can be neglected in the wide band limit. In the 
causal representation only 5 diagrams remain, in contrast to 
formulations using "naive" super-operators^ where 20 terms 
remain which partially cancel out. 



the E diagram can be performed by closing the integra- 
tion contour in the upper/lower half-plane of the complex 
plane for 77 = =p and neglecting small corrections of order 
T/D < 1: 



E(z) = / dujd 



j — rjuj — 77/U1 2 

(71) 



where Ti = T r for multi-index 1 = r],cr,r,u>. The self- 
energy E has a clear physical meaning: it is the self- 
energy one obtains in the infinite temperature limit. This 
follows from its definition since all contractions 7 = 
for T = 00, and E is the only diagram in wide-band 
limit without this contraction function. As discussed 
in Sec. B 2|) the retarded correlation function of the 
reservoirs 7 contains only spectral information (see the 
discussion of Eq. (p33")) ) and is therefore independent of 
T. Therefore E = E at T = 00. Since at T = 00 no 
energy scale matters any more, E is independent of the 
QD frequency z or any QD energy scale in the problem 
as well as the cut-off D. The dependence on T remains, 
however: in the high-temperature limit all quantum dot 



states are equally accessible by tunneling processes. This 
is described by the self-energy E. 

The action of E is very different from that of the Li- 
ouvillian L of the isolated QD: it is not super-Hcrmitian, 
as L is, but rather anti-Hermitian, 



Et 



(72) 



which simply follows from the Hcrmitian conjugation 
property of the causal representation Eq. (j5"7| . For 
T = 00 the effective Liouvillian thus reduces to L(z) — 
L + E(z) with stationary state 



P 



£1. 



(73) 



which is the maximal entropy state. This follows directly 
from the causal representation of E, Eq. (fTTj) . which 
shows that E and G share the same right eigenvectors, 
combined with the probability conservation Eq. (|58[) . By 
the same argument it is clear that in this limit the cur- 
rent vanishes: anticipating the result Eq. (|185|) from 
Sec. (|IIEp . we find for the self-energy required for the 
current equals 



E r (z) = - z ^G 1 G I | ri=r 



(74) 



where in the sum over 1 we exclude the reservoir index 
r. This also vanishes by Eq. (|55|) and therefore (I r ) = 0. 

d. Cut-off dependence and complete basis Another 
advantage of the causal representation is that the cut-off 
dependence of integrals in the individual diagrams that 
do contribute in the wide-band limit can be analyzed on 
the level of super-operators. These self-energy contribu- 
tions seem to depend on the energy integral cut-off D. 
However, using the causal structure of the perturbation 
theory one can explicitly see that such dependence can- 
cels out due to the superoperator structure of the vertices 
(i.e. the matrices multiplying the Z?-dcpcndcnt contribu- 
tions integrals vanish). The condition for this is that one 
keeps the complete basis of many body eigenstates of the 
dot Hamiltonian. To see this, however, one needs to con- 
sider the entire Liouville space, including all off-diagonal 
density operator elements and not restrict the analysis 
to only diagonal density matrix elements based on sym- 
metry properties as is often done. The idea is best illus- 
trated by considering the 1 loop contributions in Fig. [T] 
By Eq. (|71|) the 1-loop diagram with a 7 contraction is 
explicitly independent of D. For the 7 contraction of two 
G type vertices 



dujG\ 



7(770;) 



z — L — rjU) — 77/J1 



Gi- 



(75) 



the D dependence enters through the most divergent part 
of the integral, obtained by neglecting Lo in the denom- 
inator. This part is thus proportional to the superoper- 
ator 



G\Gi — — GjGi — G2G2 — 0, 



(76) 
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which is identically zero due to anti-commutation rela- 
tions (|53p (we renamed the dummy summation indices 
1 = ?/, a, r, w; in fact, only the summation over r/i is rel- 
evant here). Importantly, this argument breaks down 
as soon as many-body states have been excluded from 
the Hilbert space basis, e.g., based on their large en- 
ergy. For example, in the limit U — > oo one can exclude 
the doubly occupied QD state |2) and thereby eliminate 
Liouvillc-space elements with eigenvalues exceeding D. 
This simplifies the calculations, but for vertex operators 
projected onto this subspace the relation Eq. (f76f docs 
not hold anymore (since it is a non- linear relation). As 
a result the explicit D-dependence remains and does not 
cancel out from such expressions and the cut-off should 
be set to D ~ U if one makes this approximation. 

The above analysis can be extended to higher-order di- 
agrams and one finds that also there the D dependence 
drops out. Here one uses that diagrams containing G ver- 
tices are always expressible in terms of the cut-off inde- 
pendent E skeleton (see also Sec. pi B 4[0 . This analysis 
confirms observations made after explicit calculations of 
the kernel, e.g.— but on a very general level: the cut-off 
dependence can be completely assessed based on supcrop- 
erator algebra, without the need for explicit calculation 
of the matrix-elements of £. 



4- Finite temperature perturbation theory - elimination of 
infinite temperature contractions 

As already mentioned at the end of Sec. pi B 3|) . the 
causal structure of the perturbation theory naturally sug- 
gests to proceed in two stages. We first eliminate all di- 
agrams containing the retarded contraction 7, i.e. the 
skeleton diagram Eq. (|7ip . Since this contraction always 
occurs isolated and inside other Keldysh contractions (c.f. 
Sec. pIB 3jl ). we can use (|71[) as a skeleton diagram and 
on each propagator line resum the series, 

1 ±(t 1 y. , ,77, 

z-L-X ^\ z-L-XJ z-L~X ' 

thereby renormalizing the dot Liouvillian to 

Z = L + E. (78) 

This is illustrated in Fig.[2j As discussed in Sec. piB3cl) . 
E equals the QD self-energy in the limit of infinite tem- 
perature in all reservoirs. As discussed in Sec. pi B 21) 
the retarded correlation function of the reservoirs con- 
tains only spectral information through the energies (see 
the discussion of Eq. (|63|)). The lack of energy depen- 
dence of £ reflects that in the high-temperature limit, 
all QD states are equally accessible via transitions in- 
duced by the electrodes. Physically, one expects that the 
Keldysh contractions 7, describing the non-trivial, tem- 
perature dependent part of the distribution function^, 
should drop out in this limit. 




1 



z-ui - (£ + £) 

FIG. 2: Example resummation of diagrams with a single 7 
loop (solid blue curve) and n = 0, 1, 2, . . . skeleton diagrams 
E, resulting in 1 loop diagram with a renormalized propaga- 
tor. Since in the wide-band limit 7 contractions (dashed blue 
curve) cannot contain other vertices this can be extended to 
any irreducible diagrams with any number of 7 loops. 

In a second stage we calculate the finite temperature 
effects which arc all incorporated through the self-energy 
E(z). Its diagrammatic perturbation scries has the same 
structure as for E, but is expressed entirely in terms of 
the contraction 7, the vertex G and the Liouvillian L: 

E(z) = (79) 




Summing this renormalized perturbation theory, one ob- 
tains of course the same exact result Eq. ([55]) for the 
effective Liouvillian, which, however, is now decomposed 
in a new way: 

L(z) = L + t(z) (80) 

Since all contractions 7 have been eliminated simultane- 
ously in the above two stage transformation of the pertur- 
bation series we previously referred to it as the discrete 
step of a renormalization group procedure*^. In the An- 
derson model in the wide-band limit only one skeleton 
diagram contributes to L through E and no renormaliza- 
tion of vertex G is required to eliminate 7^ This dis- 
crete first step is a necessary preparation for the second 
step RG in which we will integrate out the 7 contrac- 
tions as well in a continuous RG-flow. In contrast to 
the G, the vertices G do not share the left zero eigen- 
vector (Zl\ with the effective Liouvillian, cf. Eq. ([86]) 
below. If they were still present during a process of con- 
tinuous renormalization, this would lead to divergences 
whenever the zero eigenvalue of L would appear in the 
resolvent (z — A — L) _1 , see also Ref . H^. The vertices G 
must therefore be integrated out before any continuous 
RG can be formulated. This is a characteristic feature 
of an RG for dissipativc systems that exhibit stationary 
states (zero eigenvectors). The causal structure of the 
renormalized perturbation theory (|79p makes clear that 
this elimination has been achieved: the zero eigenvalue 
term of L always drops out due to the presence of G's 
adjacent to the propagator (z — X — Z)" 1 . 

At this point we can already make an observation 
which is important for the construction of an explicit ex- 
pansion for the effective Liouvillian L(z) in Sec. pi C 5[) : 
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no continuous RG scheme, which reorganizes the pertur- 
bation series Eq. (|79p involving only G vertices, can ever 
generate terms of the form G...G, i.e., with a G super- 
operator standing on the far right. In the notation of 
Sec. (jll C 5|) this implies that the coefficient of the term 
\Zr)(Zr\ in the expansion Eq. (|137|) is not changed dur- 
ing a RG flow (or any non-perturbative approximation 
to the series ([75)) ). 

What was done so far can be understood as a formal 
expansion around the infinite-temperature limit as a ref- 
erence point. We start from the exact (in the wide-band 
limit) solution for the infinite temperature and then re- 
formulate perturbation theory for the finite temperature 
(including also the zero-temperature limit which we are 
going to consider) such that it explicitly contains the 
infinite-temperature limit. Note that this is not to be 
confused with the expansion in powers of T: we formally 
expand in the Keldysh distribution functions functions 
(contractions) 7. Actually one meets this expansion al- 
ready in most standard applications of the generalized 
master equation when calculating the kernel ("rates") 
such that they correspond to Fermi's Golden Rule. In 
this case one calculates only the corrections to E to the 
first order in 7. Here we reformulated the exact per- 
turbation theory around this reference point in order 
to perturbatively calculate the important higher order 
corrections to that result. To this end we performed a 
second natural step by exactly incorporating the "triv- 
ial" infinite-temperature fluctuations into a redefinition 
of the Liouvillian. We will develop an RT-RG approach 
in order to incorporate "non-trivial" finite temperature 
fluctuations also non-pcrturbativcly. This second step is 
the crucial starting point for it since this step prevents a 
serious technical problem related to the zero-eigenvalue 
from appearing, see Sec IIII Al and Sec lII C 5l 



C. Basis of the Anderson-model 
in Liouville Fock-space 

To obtain explicit RG equations that can be efficiently 
solved numerically we need to construct a basis that ex- 
ploits the advantages of (i) selection rules for the self- 
energy superoperator induced by the global symmetries 
(ii) the causal structure of E(z) and (iii) fermion-parity 
superselection rule. Again, these can be useful in per- 
turbation theory as well. 



1. Liouville space bra-ket formalism 

By Eq. ||B5J) and {73$ the QD self-energy E (E) is a 
functional of the QD Liouvillian L (L) and the vertex 
superoperators G q (G) which linearly act on the QD 
Hilbert space of many-body states. To explicitly cal- 
culate these it is convenient to introduce a "bra-ket" for- 
malism analogous to that of standard quantum mechan- 
ics: we write a 16-component supcrvector (representing 



an operator acting on the 4-dimensional Hilbert space) 
as a rounded superket A = \A) and introduce its dual 
supervector (A\» = Tr ulJ») as a linear functional act- 
ing on operators. An operator A is orthogonal to B if 
(A\B) = Tr(At.B) = (see for example Ref. |H). The 
dual vector may be written as the Hcrmitian conjugate 
of a superket: 

(A\ = L4)t, \A) = (A|t (81) 

where f is not to be confused with the supervector corre- 
sponding to the Hcrmitian conjugate of the operator A, 
A' = \A'). An operator basis \Ai), i = 1, . . . , 16 of mutu- 
ally orthonormal operators Tr(A\Aj) = Sij is complete 
if any operator B can be expanded as B — 'Yl li {Ai\B)Ai 

with coefficients (Ai\B) = Tr(A\B). Any superopera- 
tor S acting on such operators can then be expressed in 
general as a sum of 256 terms: 

16 

S=J2( A i\S\A j )\A i )(A j \. (82) 

2. Liouville Fock-space basis 

To maximally reduce the number of terms in the ex- 
pansions of superoperators (|82p we now exploit the close 
analogy of Liouville field superoperators in the causal 
representation to the usual fermionic field operators. In 
this section we first construct a suitable orthonormal ba- 
sis of operators (supervectors) in which any QD opera- 
tor can be expanded. In Sec. pi C 3[) we analyse their 
transformation properties under the symmetry transfor- 
mations of the Anderson model. In this basis we can 
then easily construct superoperator expansions compati- 
ble with these symmetries in Sec. C 5|) . 

We start from the key property (|55|) of the vertex op- 
erators in the causal representation. In bra-ket notation 
the trace operation 

Tr D . =Tr D U = 2{Z L \ (83) 

corresponds to the action of the dual on the normalized 
supervector 

\Z L ) = §1. (84) 

Therefore, by Eq. (|55]). (Zr,\ is a left zero super- 
eigenvector of G and by ([57| it follows that \Zl) must be 
a right super-eigenvector of G: 

(Z L \G 1 =0 G 1 \Z L )=0 (85) 

We can formally consider the state annihilated by oper- 
ators G as a "super- vacuum" . (Note that the "super- 
vacuum" state is the most symmetric dot operator sim- 
ilar to standard field theories where the vacuum state 
usually is the most symmetric one). This vacuum super- 
vector is proportional to the physical infinite temperature 
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density operator (cf. Eq. ([73])) with maximal von Neu- 
mann entropy S = —phi(p). It is then also natural to 
construct the corresponding right zero super-eigenvector 
of Gi whose dual is the left zero super-eigenvector of G\ : 

Gx\Z R ) = Q (Z r \G 1 = (86) 

Since Gi is a non Hcrmitian supcropcrator these eigen- 
supervectors are not simply related by Hermitian conju- 
gation in Liouville space, \Z R ) ^ {Zj^ (cf. Eq. (|8T|) ). 

The operator Z R will be found to play a key role 
throughout this work. In App. [E] we discuss its many 
interesting properties, including its similarity to Grass- 
mann numbers and its relation to the spin- and charge 
rotations of Sec. (jll C 3|) . To construct \Z R ) we now ex- 
ploit the close analogy to regular field operators: G\ is a 
creation operator in Liouville Fock-space, since its Her- 
mitian super-conjugate G\ annihilates the super-vacuum 
state \Zl) by Eq. (|85|) . The state annihilated by Gi is 
therefore simply the maximally occupied state in the QD 
Liouville Fock-spacc, starting from the vacuum^ 

i^)=n(n^J ^ ( 8? ) 

As for the regular fermionic field operators, action of cre- 
ation operators preserves the normalization: this follows 
directly using Eq. (|81[) and the anticommutation relations 
([53]): 

(Z R \Z R ) = (Z L \l[(l[G fja G r><7 ) \Z L ) 
f V n / 
= (Z L \Z L ) = 1 (88) 

Note that by Eq. ([53]) reordering the G's in the definition 
of Z R only amounts to an unimportant redefinition of the 
overall sign. Using Eq. ([5T]) we can write for the explicit 
action on any operator (denoted by •) with even fermion- 
parity 

G +CT G_ ff . = | (K, •]+ - 1 + 4 • d a - d a • 4) . (89) 
Inserting this into Eq. ([57]) we obtain 

\Z R ) = ll[(2n a -l). (90) 

Since (Z^\ is a left eigenvector of G (cf. Eq. ([55]) ) we 
immediately sec that the two zero eigen supervectors are 
orthonormal, (Zl\Zr) = 022- Therefore, it is natural to 
include \Z R ) into the orthonormal Liouville space basis. 
By successively acting on the super- vacuum state \Zl) 
we can generate more normalized, orthogonal operators. 
There arc in total 8 bosonic operators: 



Zl) 


= _ 


(91) 


\X*) 


= G+aG-^Zjf), 


(92) 


\T n ) 


= tiG^G^IZl), 


(93) 




= G+ctG- S \Zl), 


(94) 


Zn) 


= G + ^G--fG + iG~i\Z L ), 


(95) 



and there are 8 fermionic operators: 
\a+ ta ) = G +a \Z L ), 

|a+ CT ) = aG- S \Z L ), (96) 
l a +,o-) = G +a G +s G- S \Z L ), 
|a" CT ) = o-G- S G-vG +a \Z L ). 

The labeling of these basis super vectors is motivated 
by their explicit expressions in terms of the field opera- 
tors d CT ,4 ( see Eq. (|103[) - ([113[) ) and their behavior un- 
der symmetry transformations which will be discussed in 
the next section. We see that the index v of the op- 
erators \o^ a ) only has a meaning in Liouville space: it 
distinguishes v = + states with one excess excitation rel- 
ative to the super- vacuum \Zl) (super-particles), from 
v — — states with one deficit particle with respect to 
the maximally occupied super-state \Zr) (super-holes). 
By construction these 16 operators form a complete or- 
thonormal basis of the QD Liouville space. This ba- 
sis includes only one operator which has non-zero trace, 
namely \Zl) = ^1. Since all other supervectors are or- 
thogonal to \Zl), their corresponding operators are trace- 
less according to Eq. ([8"3"]) . 

We emphasize that the choice of basis supervectors 
\Zl) and \Zr) relies on a general physical property of 
the problem, the probability conservation (cf. Eq. ([6"9"]) 
and Eq. ([85]) ). The choice of the signs of the remaining 
operators, however, is motivated by considering the sym- 
metry transformations of the Anderson model. These 
will be discussed in the next section. We note that sim- 
ilar Fock-space and superoperators were introduced by 
Prosen2i see App. ([C]) . 

Before we proceed, we emphasize the necessity of work- 
ing with a complete basis for the QD Liouville space, 
which includes operators that are non- diagonal in both 
spin and/or charge quantum numbers. In perturba- 
tion theory, one can disregard all matrix elements of 
the self-energy S(z) involving the latter operators (cf. 
Sec. piB 3[) ). due to the conservation laws (|101l) . How- 
ever, the inner part of a diagram contributing to S(z) (i.e. 
a virtual intermediate state) the states are less restricted 
by the conservation laws, requiring the matrix elements 
of the vertices and the QD Liouvillian L between all of 
the above off-diagonal operators^ In the RT-RG that we 
set up below such matrix elements involving odd fermion 
parity operators cannot be avoided for S(z) as well since 
one needs to describe the renormalization of all virtual 
intermediate states, both fermionic and as well as bosonic 
ones. 



3. Irreducible transformation 
under symmetry operations 

For the total system of QD and reservoirs both the 
charge and spin component along the magnetic field are 
conserved, cf. Eq. ([2"5]) - ([2^]) . The particle number conser- 
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vation on the dot is conveniently expressed using charge- 
polarization operator 



T, = i(n-1) 



(97) 



which measures the difference of the occupation proba- 
bilities of the empty and double occupied QD states |0) 
and 1 2). It follows from Eq. and (HU) that the dot 
supcropcrators for charge polarization (cf. Eq. (1431) ) and 
spin, 



(98) 
(99) 



respectively, commute with each other as well as with the 
Liouvillian L and the self-energy S(z). Therefore the ef- 
fective Liouvillian L(z) = L + E(z) which determines the 
time-evolution of the reduced density operator conserves 
these super-observables (see App. IFlfor a derivation): 



\L(z),L T ^ 



0. 



[L(z),L s >]_ =0. 



(100) 
(101) 



Thus L{z) can be simultaneously block diagonalized with 
the superoperators for charge- and spin polarization and 
they have common super-eigenstates. In the basis of 
Sec. (|II C 3P the initial dot Hamiltonian operator Eq. ([T]) 
has the form: 

H = (26 + |)|Z L ) + ^\Z R ) + 2 + |) \T Z ) + B\S Z ) 

(102) 

In two cases the total system has a higher symmetry 
group: for e = — U/2 and /i£ = fin it exhibits full charge- 
rotation symmetry and for B — full spin-rotation sym- 
metry. These full symmetry groups are obtained by 
adding S x and S y operators to S z (cf. Eq. (|106|) ) to 
obtain the SU(2) spin-algebra and adding T x and T y to 
T z (cf. Eq. pTU|) ) to obtain another SU(2) algebra gen- 
erating "charge rotations" (see below). The construction 
of superoperators that transform in the simplest possi- 
ble way under these symmetry operations is greatly sim- 
plified by first constructing basis operators that can be 
classified with respect to the corresponding irreducible 
representations. Importantly, the Liouville-Fock basis 
(f9"Tj) - (|96p that we constructed using the causal field su- 
peroperators is already very close to a symmetry-adapted 
basis and we now merely complete the classification. We 
first group these operators according to their even or odd 
fermion parity, and refer to these bosonic resp. fermionic 
operators. We subsequently classify them according to 
the transformation behavior under the two SU(2) rota- 
tion groups as an irreducible tensor (ITO). This allows us 
to identify which supervectors correspond to diagonal op- 
erators acting on the Hilbert space, either with respect 
to the QD charge and/or spin: in perturbation theory 
one only needs the self-energies connecting such diagonal 
components of the density matrix. In the RG we set up 



below this is no longer true but it is important to single 
out this block of matrix elements of S(z), S(z) and L(z). 

We have 8 bosonic operators with integer charge and 
spin ITO ranks: 

• The zero-eigenvectors of the vertex supcropcrators 
are scalars (rank spin- and charge-ITO) with re- 
spect to both spin and charge rotations since the 
zero-eigenvalue equations Eq. ([85)) - ([86]) arc invari- 
ant under these transformations. They are there- 
fore charge- and spin-diagonal operators. This is 
also clear from their explicit form in terms of the 
Casimir operators of the spin- (S 2 = Sf) and 
charge- rotation (T 2 = £\ T 2 ) SU(2) Lie algebras: 



\Zl) 
\Zr) 



|1 

2n T n; 



(T 2 + S 2 ) 
(T 2 - S 2 ) 



(103) 
(104) 



• The generators of rotations in spin space: 

\ s *) = ^Y, a \^ = -Jill™* 

\S (T ) = S tT d a , a=t,i 



(105) 



Operators: — \S+), \Sq), \S-) transform as the com- 
ponents of a rank 1 spin-ITO (vector), i.e., simpler 
than the more familiar Cartesian components of the 
spin-operator, 

\S*, V ) = QS+) ± \3-)) , \S Z ) = ^\S ) (106) 

which satisfy [\Si), \Sj)] = ieijk\Sk)- They trans- 
form as a rank ITO (scalar) with respect to 
charge-rotations and are therefore charge-diagonal. 

• The generators of rotations in charge space: 

l T ") = 7fEl^) = 72("- 1 ) ( 107 ) 



|T+) = 44; 
\T-) = 



(108) 
(109) 



Operators: — |T+), \Tq), |T_) transform as a rank 
one ITO (vector) under rotations in charge space, 
and MS M rank ITO (scalar) with respect to spin 
rotations. Therefore they are all spin-diagonal. 
These are more convenient than the Cartesian com- 
ponents 

\T*,y) = 2®(\T+)±\T-)), \T*) = ^\To) (HO) 



V2 1 



satisfying the SU(2) algebra [|T<), \T 3 )] = iey fc |T fc ). 
In contrast to the spin-operators, here the indices 
x, y, z are not related to the axes in the real space, 
but merely label the components of the SU(2) gen- 
erators. 
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In the bosonic sector we can thus use either 

\x<r) = n„-$l = j s [\T ) + <T\S )] (111) 

or | So) and |Tq) as basis vectors. Only the latter two are 
adapted to charge- and spin- rotation symmetry, but the 
former two allow for greater notational simplicity. Both 
will be used. 

The other half of the QD Liouville space is spanned by 
8 fermionic operators with half-integer charge and spin 
ITO rank: 

• The fermionic operators, denoted by act in 

three subspaces of dimension two labeled by a = ± 
(spin), 77 = ± (particle- hole) and, additionally, 
v = ±. These basis operators are ITOs of rank i 
with respect to both charge (77) and spin rotations 
(a). For v — + these arc the rank ^ spin-ITOs con- 
structed from creation and annihilation operators: 

K a ) = ^4, i<j = ^d s , (112) 

\a+ <a ) = V2\Z R a+ a ), |a" CT ) = V2\Z R a± a ). (113) 

All these operators are both charge- and spin off- 
diagonal. We give their above explicit form to em- 
phasize that there is an additional set of fermionic 
operators in the charge off-diagonal subspacc which 
is linearly independent of the standard field oper- 
ators d a and S a , see the discussion of the index v 
in Eq. (|96|) . Noting the property (2Z R ) = 1 we 
find the following explicit relation between these 
two sets: 

= 2Z R \al a ) (114) 
where as usual u = — v. 

The above simple transformation behavior of the basis 
supervectors motivates all the relative signs we antici- 
pated in writing Eq. (|?T|) - (|^5)) . The basis is therefore 
completely fixed up to irrelevant phases by general phys- 
ical properties (probability conservation, symmetries) to- 
gether with the supervector normalization. Table U sum- 
marizes the transformation properties. 

Finally, we note the transformation behavior under the 
action of the Hcrmitian conjugation K (see Eq. (1701) ). For 
all diagonal basis operators D = Zl, Z Rl Xcr, So, To 

K\D) = \D) (115) 

whereas for the off-diagonal operators 

K\S ± ) = |S T ), (116) 
K\T ± ) = \T T ) (117) 
K\al a ) = -vw\a\ s ) (118) 

4- Expansion of the vertices 

The calculation of bra-ket representation of the vertex 
operators in the basis Eq. ([9~Tl) - (l9"6"]) reduces entirely to 



Operator 


fermion 




I -1 ±KJ 

(rZ)Ti\( lUflPYl 
lICLllix, 1J1UCA 1 






\Zi) 


+ 


(0,0) 


(0,0) 


i = 


L,R 


\ S m ) 


+ 


(l,m) 


(0,0) 


m — 


0,±1 


\T m ) 


+ 


(0,0) 


(l,m) 


m = 


0,±1 






<- -) 


(- 2 ) 


77,(7 


= ±1 



TABLE I: Fermion-parity and irreducible transformation be- 
havior of basis operators under spin and charge rotations. 
Schematically denoting these operators by \s, m s ; t, mt), these 
transform as (i) spin-irreducible tensor operators (S-ITOs) 
with rank s and index m s , L Sz \s, m s ;t, m f ) = m s \s, m 3 ;t, mt) 
and L s± \s,m s ;t,mt) = \Js(s + 1) — m s (m s ± l)|s, m s ± 
l;t,m t ). and (ii) as charge-ITOs (T-ITOs) with rank 
t and mt, L Tz \s, m s ; t, mt) = m s \s,m s ;t,m t ) and 
L T± \s,m s ;t,m t ) = y/t(T+ 1) — m t (m t ± l)\s,m s ;t,m t ± 1). 
All zero-index S-ITOs (T-ITOs) correspond to operators act- 
ing on Hilbert space that are diagonal in spin (charge). 

the calculation of 

G +a =\a+, a )(Z L \+o-\Z R )(aZ, s \ (119) 
+ a\T + )(a+^\ + |a" CT )(T_| + |«- a ){ Xs \ 

- <t|x«0«*| + *\S ff ){a± tir \ - \oT^ s ){S s \ 

This is easily performed using the 2nd quantization tech- 
nique that we introduced in Liouville-Fock space, i.e., 
using algebra rather than the explicit matrix represen- 
tations of Sec. (jll C 3[) . All other vertices follow from 
general relations. First, using the transformation under 
Hermitian-conjugation (see Eq. (1701) ) 

<V - (-l) Ln+l KG na K (120) 

the result for opposite charge index 77 = — follows using 
the transformation properties Eq. (|115l) - (1118[) : 

G-a = - o\at. s ){Z L \ - \Z R )(a-J (121) 

- a\a^)(T + \ - |r_)(a+J + a\aZ, ff )( X *\ 
-\ Xa )(a+J-\S s )(a+J+cr\aZ, a )(S a \ 

The vertex supcroperators G^ CT can be obtained from 
Gf/,<T by Hermitian conjugation in the Liouville space us- 
ing Eq. (157)) . In the Liouville bra-ket formalism this sim- 
ply means that we can interchange bra and ket vectors in 
the expansions Eq. (|119|) and Eq. (|121[) to obtain G Ta . 

We point out that under the 2-loop RG-flow to be dis- 
cussed in Sec. (jIII[) the structure of the vertex operators 
is changed and the above relations cease to hold (G will 
be modified whereas G is not changed, implying that 
Eq. (j57|) breaks down). However, in Sec. (jIV C| we show 
how a such vertex corrections can be incorporated effec- 
tively into the flow of the effective Liouvillian only, al- 
lowing us to work with the "bare" vertices with the nice 
properties discussed above. One property of the bare 
vertices that remains valid under the RG-flow is 

(o+|Gi oc (Z L \ or 0, (122) 
G\a~ ) cx \Z R ) or 0, (123) 
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Eq. (|122[) follows since \a+ a ) is obtained from the vacuum 
by the action of a single creation superoperator (cf. 
Eq. Using Eq. (HZJ) 



In contrast, 



(124) 



When inserted in the left hand side of Eq. (|122l) this 
G can be anticommuted past the G using Eq. ([5"3"]) and 
the super vacuum property (Zl\G = (see Eq. ([85]) ). 
Analogously Eq. (|123l) follows by noting that Eq. ([91)]) 
for the la^g-) can be rewritten as a single destruction 
superoperator acting on the maximally occupied state in 
Liouvillc-Fock space: 



\a- M ) = a^)/i GrhM \Z R ) 



(125) 



Commuting the G to the left and using G\Z R ) — (see 
Eq. ([55)) ) on the right of each term we obtain Eq. (|123[) . 



5. Expansion of the effective Liouvillian 

a. Causal structure By the general properties (|85|) 
and (j86p the vertices must have an expansion of the form 
(confirmed by the explicit results Eq. (fTT9l) - ([T2T[) ): 



G=..\Z R )(.\ + ..\.)(Z L \ + 
G=..\.)(Z R \ + ..\Z L )(*\ + 



(126) 
(127) 



where the remaining terms involve neither Zl or Z R . 
Therefore the terms in the expansion of the effective Li- 
ouvillian involving these vectors are strongly restricted. 
Combined with the general causal structure of the per- 
turbative series, i.e., the way G and G can appear, this 
imposes further -constraints (cf. Sec. ([II B 3| ): 

• Terms of the form |Zl)(»| are prohibited by prob- 
ability conservation, since otherwise the trace con- 
dition (1691) would be violated: 



(Z L \L =(Z L \L(z) 
= {Z L \t 



(Zl\Z(z) 
(Z L |E(z) = 



(128) 
(129) 



This is guaranteed by the causal structure which 
requires that the leftmost vertex is always of the 
type G, with expansion Eq. (|126|) . 

• Terms of form |»)(.Zr| can only appear due to the 
diagrams collected in E = — ijTiG\Gi. Expand- 
ing Eq. ([71]) in the basis ([II])- ([M]) one finds that 
only the term \Z R )(Z R \ with coefficient — i4r can 
appear, see also Eq. (|135[) . Importantly, this im- 
plies that \Z R ) is a right eigenvector of both L as 
well as the exact effective Liouvillian L(z) and the 
kernel E(z) , see discussion of Eq. (|137[) : 

L\Z R )=L{z)\Z R ) = V{z)\Z R ) (130) 
= t\Z R ) = —4:iT\Z R ) (131) 



2(z)\Z R ) = 



(132) 



Note that this eigenvector and eigenvalue are inde- 
pendent of the QD frequency z. 

The term \Z r )(Zl \ is not forbidden by general con- 
siderations. However, such terms always drop out 
in the calculations of the transport current that 
interest us here (see Sec. piEjl ). This happens be- 
cause in all required expressions the renormalized 
Liouvillian L, parametrized as (j!37j) . is evaluated 
between two G vertices cf. Eq. ([79"]) . Therefore, 
by Eq. ([T26]) the term \Z R )(Z L \ with coefficient ( 
always drops out. We emphasize, however, that £ 
does enter into the stationary state (cf. Eq. (|152l) ) 
and other physical quantities, such as the average 
dot energy, and may therefore be important for, 
e.g., thermal transport problems. 

Terms of the form |a+ CT )(»| and | •) (ck~ CT | can ap- 
pear only in the bare Liouvillian L = [H, •] or the 
infinite-temperature kernel E = — i^FiGiGj, but 
not in the non-trivial kernel E(z): 



(a+Mz) = 
E(z)K„) = 



(133) 
(134) 



Both relations follow from the fact that £(z) con- 
tains only vertices of the type G in expansion 
Eq. ([7J|]). Eq. (fT3"3"]) follows from Eq. (p^j) applied 
to Eq. ([79]) and using (Z L \L = and (Z L \G = 0, 
the vacuum property (see Eq. ([85]) ). Analogously 
Eq. (|134|) follows from Eq. (|123|) using our general 
result L\Z R ) — —4iT\Z R ) (see above) and then 
G\Z R ) = (see Eq. ([56])). Eq. ([133]) -([134]) allow 
us to make general predictions about the excita- 
tion spectrum of the exact dot Liouvillian, £(z), 
see Sec. (lITDlfl) and (IIIIB4I) . 

b. Spin and charge rotation symmetry We now first 
expand the infinite-temperature self-energy E and the 
renormalized Liouvillian L in the basis ([9~T]) - (j96[) . Sub- 
stituting the above bra-ket expansions of the supcropcr- 
ators Gi and Gj into Eq. (|71[) we get: 

£ = -iT[4\Z R )(Z R \ + 2^2\ X< r)(x<r\ 

cr = ± 

+ 2^|T CT )(T a |+2^|5 ff )(5 ff | 



cr = ± 



a— ± r\— ± v— ± 



(135) 



Clearly, E is explicitly anti-Hermitian in the superoper- 
ator sense (cf. Eq. ([72"]) ). Combining this with the bare 
dot Liouvillian, obtained by expanding the commutator 
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L = [H, •] of (fT02|) we obtain 



-^a|a=oo 



L = L + T, 



HT\Z R ){Z R \ - 2iT Xo 
J2 [{<rB - 2iT)\S <r )(S <r \ + (a(U + 2e) - 2iT)\T a ){T a 



E [v 

^"E \ a v,<r)( a v 



€+ 2j + <T f -*( 2 -^ r 



(136) 



where for reference we introduced the notation 
£a|a=oo := L of Sec. . All non-zero eigenval- 
ues of L have negative imaginary parts, thereby au- 
tomatically regularizing all resolvents that can appear 
G(z — L — X) _1 G in the perturbation theory for 
for z — > iO. This can be seen explicitly since all terms 
are already in diagonal superoperator form, with the 
exception of the odd-fermion terms (the last two lines 
in Eq. (|136[0 whose eigenvalues are given below (set 
AF~i + = in Eq. (|164l) ) . Note that the right zero eigen- 
vector of Za =00 is \ | Z L ) - it is the only basis vector miss- 
ing in Eq. (|136|) - in agreement with the result (|73|) in 
Sec. (jll B 4[) . We further note that the infinite tempera- 
ture kernel £ (|135[) contributes terms to Eq. (|136l) that 
are diagonal in the index v, whereas contributions off- 
diagonal in v arc produced in Eq. (|136j) by the Coulomb 
interaction included in the bare dot Liouvillian L. Wc 
will show in Sec. B 4[) that the continuous RT-RG to 
be set up in the next section produces only i^-off-diagonal 
contributions in the effective Liouvillian L(z) with impor- 
tant consequences. 

We can now write down the exact form of the QD 
effective Liouvillian L(z) taking into account all general 
restrictions that we have derived above. In the most gen- 
eral case that we consider only restricted spin- and charge 
rotation symmetry about an the z-axis: this implies that 
the effective Liouvillian must be a sum of superoperators 
that (i) transform as an irreducible tensor of any rank 
but with index zero with respect to both the charge- and 
spin-rotation group (i.e. by pairing only bras and kets 
of basis supervectors with the same charge and spin in- 
dices) (ii) preserve the fermion-parity (i.e. by avoiding 
combinations of fermion and boson kets and bras). Us- 
ing Table U which lists these transformation properties 
of the basis supervectors (f!H"j) - (jM)) . we can readily con- 
struct the most general form of superoperators of this 
kind that are furthermore compatible with the causal 
structure of the perturbative series (|68|) . In Table |TT] 
we have classified all these superoperators according to 
their irreducible transformation properties under the full 
symmetry group of both spin and charge rotations. This 
makes it easy to impose further restrictions on the ex- 
pansion coefficients in cases of higher symmetry (B = 
and / or e = -U/2,(i L = Mi?), sec Sec. (|IID4|) . The 
most general form of the exact QD effective Liouvillian 



Superoperator 


S-ITSO 

(rank, index) 


T-ITSO 

(rank, index) 


\Zr)(Z r \, \Z r )(Z l \ 


(0,0) 


(0,0) 


E \Sm)(S m \ 
m = 0,±l 


(0,0) 


(0,0) 


E \T m )(T m \ 

m = 0,±l 


(0,0) 


(0,0) 


\So)(Zl\, \Z R )(S a \ 


(1>0) 


(0,0) 


?7i i jS™ ) ( S™ i 

/ j ""I 1 -' 777 / \ 777, | 

m = ±l 


(1 0) 


(0 0) 


\T )(Z L \, \Z R )(T \ 


(0,0) 


(1-0) 


E rn\T m )(T m \ 
m=±l 


(0,0) 


(1,0) 


|So)(T |, \T ){S \ 


(1,0) 


(1,0) 


e ^<vk ct )k; ct | 

cr ,77= zb 




(r„,0) 


TV = 0, 1 


r v = 0, 1 


E (Sm 2 -2)\S m )(S m \ 

m=0,±l 


(2,0) 


(0,0) 


E (3m 2 -2)|T m )(T m | 

m = 0,±l 


(0,0) 


(2,0) 



TABLE II: Irreducible tensor superoperators (ITSOs) of dif- 
ferent rank but with (i) zero index with respect to both spin 
and charge rotations and (ii) satisfying the causal structure 
constraints (cf. Sec. (|II C 5 all ). The general effective Ander- 
son Liouvillian (|137p is a linear combination of all of these, 
where the coefficient of \Zr)(Zr\ is always fixed to — 4ir in 
the wide-band limit, see Sec. ()II C 5 a|) . In the special limits 
of higher symmetry only (0, 0) S-I TSO re sp. T-ITSOs can 
appear in this expansion, see Sec. (|II D 4|) . The ITSOs are 
constructed by standard angular momentum coupling. For 
this one takes the supervectors in Table |IJ denoted schemat- 
ically by \s,m a ;t,mt) where s,t and m 3 ,mt are the rank 
and index with respect to spin and charge rotations in Li- 
ouville space. Then one constructs conjugate bra supervec- 
tors which transform with the same rank and index: these 
are (_ i) s - m »+'- m * ( S) —m s ;t, —m t \. Coupling these kets and 
bras with Clebsch-Gordan coefficients gives the superopera- 
tors transforming with definite rank and index with respect 
to spin and charge rotations. 



L(z) = L + S(z) = L + S(z) then reads 
iL(z) = 

4T\Z R )(Z R \ + C\Z R )(Z L 



(137) 



e+E 



M a \T a ){T a 



E a \S a )(S a \ 



\^){Zl\ 



where here T = | EJ r L r and F, ] cr are superoperators 
acting in the 2-dimensional v space spanned by a ). 

This is a central result of the paper and before dis- 
cussing it the coefficients in detail we point out its impor- 
tance. By exploiting only general properties we have re- 
duced the number of terms from 256 (cf. Eq. (|82|)) down 
to only just 30. The key simplification comes from using 
the causal field superoperators to construct the Liouville 
Fock-space. Moreover, the general Liouvillian can be eas- 
ily diagonalizcd as we show in the next section. Further- 
more, because of its general nature, the parametrization 
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(|137p is useful in other applications than those considered 
here and may be extended to more complex Anderson- 
type models. For example, we note that \Zr) is always a 
right eigen supervector of the effective Liouvillian decay- 
ing with rate 2mT where m is the number of electrodes 
attached to the dot (for our case m = 2). This mode was 
recently also found in a study^S investigating the time- 
relaxation of the density matrix of the Anderson-model. 
It was observed that this mode, appearing in 1-loop per- 
turbation theory, is not affected by 2-loop corrections. 
Our work generalizes this result: the eigenvalue of the 
right eigen supervector \Zr) is not affected by any higher- 
order corrections. We also see how this relies on assuming 
the wide-band limit. Both insights directly rely on the 
causal representation of the field superoperators. The 
time-dependence will be discussed further elsewhere^. 

We now discuss how the expansion coefficients for L(z) 
are incorporated in Eq. (|137j) : 

• The choice of the first coefficients of the terms in- 
volving the vertex zero super-eigenvectors \Zl) and 
\Zr) is based on the general properties of the per- 
turbative series (cf. Sec. (|II B 3|l ). 

• A supervector in the two dimensional x-subspace 
spanned \xc) 



(138) 



• An independent vector in the corresponding dual 
space 



(<i>\ = £<Mx<> 



(139) 



• A supcropcrator acting on the x-subspace 

£ = Ix.Xx.' I = £ &Xi (wo) 

a,a 1=0,1,2,3 

Here, in the second line, the matrix £ CT)Cr ' was de- 
composed in the standard basis (Tj)^' of the unit 
(i = 0) and three Pauli matrices (i = 1, 2, 3), giving 
another x-subspace basis of superoperators: 



a.a' 



(141) 



• Four superoperators acting on the 2 dimensional 
a r , !(7 -subspaces spanned by |a^ j0 .) 



Frj.a — ^ t F na \a ria )(a rilT 
= V F l a 1 



(142) 
(143) 



1=0,1,2,3 

with unit and Pauli- vector superoperators 
for each fixed a = ± and r\ = ±. 



It is convenient to use both the four-vector as well as the 
2x2 matrix representations for the superoperators £ and 
J? 

All the above expansion coefficients depend on the QD 
dot frequency z (not written) and satisfy the following 
conjugation relations, deriving from the Hermicity con- 
dition dZQj): K(iL(z))K = iL(-z*). For the charge- and 
spin-diagonal operators these are 



(-Z*), 

£(*) = £*(-**)> 



^(z)=^*(-z*), 



(145) 



implying that these are real only for zero frequency z = 
iO. For the charge- or spin- non-diagonal operators we 
have 



F^'{z) = vv'F^f{-z% 
M a {z) = Ml{-z*), 
E a (z)=EU-z*). 



± 



(146) 
(147) 
(148) 



Therefore at finite dot frequency z all parameters are 
in general complex and all 2 x 2 coefficient matrices are 
non-Hcrmitian. 

Note that in Eq. (|137|) we have parametrized iL(z), 
rather than L(z), i.e., including the imaginary factor i. 
This anticipates the application to the RT-RG where 
a renormalizcd version of the Liouvillian L appears in 
the final RG equations only the combination iL (cf. 
Sec. (|IIIC ip and (|IIID ![) ). Finally, we emphasize that 
the simplifications that led up to the parametrization 
Eq. (TT37) remain valid for the RT-RG: Eq. <j69j), as well 
as Eq. (|85jl and Eq. (|86| do not change under the contin- 
uous renormalization as will be shown in Sec. (IIIII). 



D. Effective Liouvillian 



1. Spectral decomposition of L(z) and L 

Above we have reduced L(z) to block diagonal form as 
far as possible by using symmetry and general properties. 
In the Sec. (|IIIj) we will see that L(z) is closely related to a 
renormalized version of the QD Liouvillian L that we will 
denote by La- To make this clear we have to anticipate 
a result: La interpolates between L and L(z) as the flow 
parameter A varies from ootoO: L(z) = L + f dA^j^. 
This is done by redistributing diagrams of E in Eq. ([50)) 
and including a fraction of them into a redefinition of the 
Liouvillian L. By construction S(z) is thus decomposed 
into "pieces" ^j^- with the same matrix structure that 
are accumulated during the flow. At the end of the flow 
La equals L(z). Therefore La has an expansion of the 
same form Eq. (|137[) as L(z). For notational simplicity 
we denote the expansion coefficients of La by the same 
variables as for L(z). In cases where this leads to con- 
fusion the coefficients of La arc distinguished from those 

of L{z) by indicating their A dependence, e.g. vs. 



18 



F"g / , which can, however, often be omitted. All results 
obtained in this section thus apply to both the exact L(z) 
as well as the renormalized L\. The explicit form of La 
needs to be inserted in to perturbative expressions in ex- 
pansions involving resolvent superoperators to evaluate 
RG-equations. 

We therefore need to completely diagonalize L\ such 
that it can be expanded into its eigen projectors P k = 
\X k ){X k \ = (P k ) 2 : 



(149) 



where the sum runs over the labels i of the eigenvalues. 
Here (A*| and \X Z ) are the left and right eigen supervec- 
tors of La for the same eigenvalue A fc : LAP k = P La = 
\kpk Using this complete and orthogonal set of pro- 
jectors one can then evaluate resolvent superoperators in 
Eq. JTSJ) explicitly: 

.G — —G... = ... J2 I u SpiG - ( 15 °) 



z- X -L A 



X - X 1 



The diagonalization of Eq. (|137[) is also useful for higher- 
orders of (renormalized) perturbation theory Eq. (|68[) 
and f|T9[) : when expanding the QD L and La in the 
form (|137[) it is directly adapted to all symmetries of the 
problem and one can efficiently construct explicit matrix 
representation of the self-energies £(2) and E(z), respec- 
tively. 

Applying the above spectral decomposition to the con- 
tinuous RG of Sec. (jllip involves two assumptions that 
should be pointed out here. First, we always assume that 
the zero eigen value of La is non-degenerate, correspond- 
ing to the unique stationary state for the density opera- 
tor. This is always found to be the case for the numeri- 
cally calculated RG flows discussed in Sec. However, 
in principle it may also happen that two (or more) non- 
zero eigenvalues of La become degenerate during this 
flow. If this is the case, and additionally the superma- 
trix the La has nonzero elements on its diagonal in the 
normal Jordan form in the degenerate subspace, then no 
complete cigenprojcctor basis exists. For the Anderson 
model this presents no crucial complication: the eigen 
basis can in principle be circumvented for the calculation 
of the 2 dimensional superoperators. However, numeri- 
cally we never meet such a problem in the application of 
the RT-RG theory presented below. 



2. Eigenvalues, super-eigenvectors and the stationary state 

We now diagonalize L(z) or La, parametrized as in 
Eq. (|137p . explicitly in each of its block-diagonals. We 
start with the block spanned by the charge and spin di- 
agonal bosonic operators \Zl), \Zr), \Sq), \Tq) which con- 
tains the stationary state: 

• Eigenvalue X Zl 



0: 



P z - = 2\p)(Z L \ 



(151) 



with the stationary density operator 



sr 



Zr) (152) 



We note that the coefficient C appears only in the 
stationary state and the next projector P Zr but 
not in any other eigen projector or eigenvalue. 

Eigenvalue X Zr = —iAT. 



pZn = ( 



1 



'4T(4T-0 



^)\Z R )(Z L \ 



(153) 



1 



\Zr){Z r \ + ^\Z R )(Z L \ + 



Eigenvalue X x± with projectors: 



P 



3,4 



px,± 



(0IP X± I^) 

(\x± - AT) X x± 



\Z R )(Z L 



(154) 



\x,± _ 4r 



\Z R )(}\Fx* + 



-h 



\x,± 



P x± \i>)(Z 



Here the eigenvalues X x± are determined by first 
diagonalizing £ in the x-subspace, i.e., finding 



(155) 



Since ^ a ,a' is a 2 x 2 non-Hcrmitian matrix, this can 
be expressed in the complex vector £ = {£,1,^,2,^3) 
and complex coefficient £0 (cf. Eq. (11401) ): 



(156) 
(157) 



X x ' a 



2v^ 



The remaining block-diagonals of L acting on bosonic 
subspaces are 1-dimensional: for a = ± 



X T ° 



-iE„, 



P T ° = \T a ){T a \, (158) 
P s ° = \S a )(S a \. (159) 



Here the square root of the complex argument is defined 
such that the branch cut lies in the lower-half complex 
plane since effective Liouvillian must be regular in the 
upper half-plane according to Eq. (|3"2"|) . (|35|) . 

Finally, in the four two-dimensional fcrmionic sub- 
spaces, labeled by a I) (T , for fixed r\ and a, the eigen val- 
ues and projectors can be calculated as for the bosonic 
X-block: 



A"*"* = -i[F* r ± 



T/<T J 



P a i, 



—or ± — ■ — =^=— 

1 2JF 2 



(160) 
(161) 



where the coefficients F® a and F, ;CT = (F^ a , F^ a , F^ a ) are 
again complex (cf. Eq. f| 143[) ^) . 
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We note that in principle is possible that £ 2 = while 
£ 7^ . In this case the supermatrix representation of x 
has non-zero diagonal element in its normal Jordan form. 
In this case x has no eigenbasis and Eq. (|156[) does not 
apply. Still, the required matrix valued functions of x can 
be calculated using the Hamilton-Cayley theorem. The 
same remarks apply to Eq. (|160[) . However in practical 
applications we never meet such a situation. We also note 
that during the RG flows discussed in Sec. (|III|) we never 
run into points where A x,± = —i4T and one therefore 
need not worry about vanishing of the denominators in 
Eq. (|154|) or the existence of the inverse of (4E — £) in 
Eq. (TT531 . 

An important simplification applies to the first four 
bosonic projectors P Zl , P Zr , P 3,4 that involve the 
vectors \Z L ) and \Z R ), analogous to the corresponding 
terms in expansion (|137l) of the effective Liouvillian (cf. 
Sec. pi C 5 a|) ). When inserting projectors into Eq. (| 150[) 
between two vertices G 

• the projectors P L , P Zr give no contributions 

• the projectors P 3,4 only contribute through the first 
term P x ^ in Eq. (fT54| . 

As a result in all applications below we can replace 
Eq. (|149|) with the simpler expansion 

L A ^ (162) 
X x P x + ^ \\ S *P S « + \ T "P T - + Y A Q "'"P Q "-' 

a rj 

Here we leave implicit the sum over the two eigenvalues 
A x,± and A Ql, ' CT,± in the x an d a^cr-subspace, respectively. 
A crucial stability requirement for the RG in Sec. (|III[> is 
thereby explicitly satisfied: the zero eigenprojector (|151l) . 
corresponding to the physical stationary state, never ap- 
pear in the resolvents. 



Fermionic excitations - spectral decomposition of £ 



The expansion (|1 37[) can of course also be applied to 
£(z) = L(z) — L as well. This however, involves ad- 
ditional simplifications causing certain terms appearing 
in the expansion of L(z) (and La) to be absent. First, 
due to Eq. (|132[) the term \Z R )(Z R \ is missing. Sec- 
ondly, due Eq. (|133l) - (|134[) deriving from the causal struc- 
ture (cf. Sec. pi C 5 a| ) most of the coefficients of the 
fermionic sector of £ , denoted by — iAF^ , vanish: 
AF^£ = S u -S u ' :+ AF~ a : + . We can express the coeffi- 
cients matrices F^ of L{z) = Z + £(z) using Eq. (|136l) . 



m terms of AF va ! + : 



f F +,+ F + -\ 

\ r}<7 r/cr / 

_f v{e+ U )+a B_ iT 



(163) 



u 

2 



U 
2 



iAF: 



Converting to spherical coefficients F® a and F na and us- 
ing Eq. (|160p we find 

X^'"'* =vU+ — \ +a—-2iT (164) 



, U 2 ra .UAF-^ + 

± v\ r 2 - % — 

'V 4 2 



v 



LLfy 1 _|_ 7 To' 3 -iAF"- + »" 

rj,a _j_ 2 r),cr ' LL lx r\.a "'- XJ - rj,cr ^rj.cr 



(J2 



T 2 



■ uaf: 



+ 



(165) 



We thus find that the functional form of the fermionic 
eigenvalues is severely limited by the from the casual 
structure of E: all the frequency dependent renormal- 
ization effects enter into the projectors and eigenvalues 
solely through the four complex numbers AF~^~ in the 
four a V a blocks (r),cr = ±). 

For an isolated QD decoupled from the reservoirs, V = 
0, we have t = = £(z) = and thus AF~'+ = 0: 

A a "'- + =r) (e + U) + aB/2 (166) 
A a "<"- = r/e + crB/2 (167) 

This is the spectrum of many-body energy excitations 
when adding a single electron, starting from either an 
empty QD (A Q ' I ''" + ), or a singly occupied QD with spin 
a (A" 11 "' - ). This is confirmed by the eigenprojectors of 
T ^ for finite U 



W,0)(a,0\ 
\2,a)(2,a\ 



(168) 



where \a,0) = \a)(0\ = (1/2 - Z R )S a and \2,a) = 
\2)(a\ = <t(1/2 + Z R )d\. These are the virtual inter- 
mediate states and energies that enter through L into 
the perturbation expansion (|68p for 

For any finite coupling T but infinite temperature, 
£ = and therefore again AF~^~ = 0. Now, however, 

£ oc T 7^ and the eigenvalues obtained from Eq. (|164[) 
depend qualitatively on the interaction strength U. For 
U < 2T the coupling to the reservoirs T adds differ- 
ent imaginary parts to the degenerate real eigenvalues 
rj(e + U/2) + aB/2 whereas for Z7 > 2T it adds a uniform 
imaginary part — iT to these eigenvalues while differently 
shifting their real parts. This dependence on U is plotted 
in Fig. [3] and shows a bifurcation at U = 2T. These are 
the energies and projectors that enter through L into the 
renormalized perturbation theory (|79p . 

Finally, for both finite coupling T and finite temper- 
ature T, the complex coefficients AF~'+ are non-trivial 
functions that need to be calculated, e.g., either perur- 
batively or using the RG, see Sec. pil B 4[) . However, 
even in this case Eq. (|164[) provides an exact relation: 
the average of the complex particle and hole excitation 
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FIG. 3: Fermionic excitation energies and widths of the 
infinite temperature Liouvillian L plotted as function of the 
tunnel coupling strength F/U for e = B = 0. The energy 
and width are given by the real and imaginary parts of the 
fermionic eigenvalues A""'"'". (Eq. (|164|) with AF~^ + = 0). 
The real parts for v = + (red) and v = — 1 (green), respec- 
tively, are given by the full lines and the imaginary parts are 
indicated by the shaded width of the level with the corre- 
sponding color. For U > 2F the excitations have different 
energies - split by U - but with the same width 2F whereas 
for 2T > U they have the same energies but different widths: 
for T > U -ImA a "'<" ± « T (red) 3F (green). 



eigenvalues in each a^-block is independent of AF V ^ 
and thereby, also independent of the frequency z: 



1 

2^ 



A' 



= V 



a 2iT 

2 



(169) 



Physically speaking, both the renormalized energies of 
single particle fermionic excitation energies (real part) as 
well as their decay rates / broadenings (imaginary parts) 
lie symmetric with respect to the above average values. 
E.g. if the particle excitation broadens, the hole excita- 
tion sharpens up and vice-versa. 

Finally, we can infer an important physical stability 
constraint on the functions AF~^~: they must be such 
that imaginary part of the square root is less than 2r 
for all A. Otherwise, inverse Laplace-transforming L{z) 
to time-space would yield terms that diverge for t — > oo, 
which is unphysical. This restricts the maximal excita- 
tion widths in the fermionic block: the negative imagi- 
nary part of eigenvalue A" cannot exceed the value 4r. 



4- High-symmetry stationary states 



As a cross check on the results of Sec. (jll D 2| . we now 
analyse the stationary density operator in the two special 
parameter regimes where the model has a higher symme- 
try than in general. In both cases the superoperator £ 
has no off-diagonal terms in the basis of \Sq) and |To), 



because there are no scalars under spin or charge- 
rotations which contain such terms as components, see 
Table ITT1 If cither symmetry is broken, charge or/and spin 
rotations allow for mixing terms |To)(S'o| and |S'o)(T'o| in 
Eq. (|137[) . For this reason we work in the general case 
with the \xa) basis (jllip . 

a. Full spin-rotation symmetry, B = 0. At zero 
magnetic field all terms in the Liouvillian must be zero 
rank ITSO (scalar) with respect to spin-rotations, i.e., 
the terms |S , )(S'o| and |S'±)(S'±| must have the same co- 
efficients, E + = E_ = £ss and the coefficients of the rank 
1 and 2 spin-ITSOs must vanish, i.e., in the x-subspace 
\<j>) = <jy T \T a ) and = 4> T \T ) with real (f> T , ip at z = iO 
(cf Eq. (fT45l) ) 

iL A = ... + T |^)(T O | + Vt|T )(Z l | (171) 
+ M + |T + )(T + |+M_|T_)(T_|+fr T |T )(T | 

+ £ss( \S+)(S+\ + \S-)(S.\ + \S )(S \) +... 



See Table \U\ where the ITSOs are listed. The expression 
for the stationary state p obtained using Eq. (|152[) is 
then independent of £ss it is therefore explicitly invariant 
under spin-rotations, as required: 

P = -^T\T ) + l\Z L )- C -^ TT \Z R ) (172) 

b. Full charge-rotation symmetry, e = —U/2 and 
= pr. At the particle-hole symmetric point, as sim- 
ilar reduction must take place: here M + = M_ = £tt, 
\<f>) = <f>s\So) and |^)=Vs|S ): 

il A = . . . + (j> s \Z R )(S \ + 4>s\S Q ){Z L \ (173) 
+ Ztt(\T+)(T + \ + |T_)(T_| + \T )(T \ 

+ E + \S + )(S + \+E_\S-)(S-\+tss\So)(S \ 

and the stationary state is explicitly invariant under 
charge-rotations: 

P = -^s\ ^o) + -| Z L ) — \Z R ) (174) 



E. Current superoperator and its irreducible 
self-energy 

Our main objective is to calculate the stationary cur- 
rent which flows through the QD. Having set up the per- 
turbation theory formalism for the density operator, the 
expression for the average current can now be compactly 
derived. Moreover, we give a general proof that in general 
at zero bias the current vanishes, as it should, indepen- 
dent of the way the self-energy is calculated. 

The current flowing out of reservoir r = L,R is 
obtained from the Heisenberg operator of the time- 
derivative 



Z = Ztt\T )(T \+£ss\S )(S \ 



(170) 



I r = rf = -i\H % 

dt L 

= -i[V r ,n r ] 



(175) 
(176) 
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Note that there is no summation over the electrode index 
r, c.f. Eq. (fTU|) . The expectation value of the current 
operator can be expressed in superoperators using the 
cyclic invariance of the total system trace: 



(r>(t)=TrTr(r P (i)) 



D R 

-zTrTr 

D 



Tr(L ir ( 

R V 



'(t-to) 



p(to)p* 



(177) 



We note that observable averages involve anticommuta- 
tors of the corresponding operator (see e.g., Ref . l4ll) : 



■[F 



(178) 



This is in contrast to time-evolution superoperators 
which involve commutators of the Hamiltonian operator. 
This difference is exploited below. If one uses Eq. (|175[) 
in Eq. (|177|) the evaluation of the reservoir trace is unnec- 
essarily complicated, since Eq. (|175[) involves two opera- 
tors acting on the reservoir (V r and n r ). Here we proceed 
differently : we first use that the tunneling through junc- 
tion r conserves the particle number of the dot and the 
reservoir r: [n r + n, V r ] = and eliminate one electrode 
operator: 



F = i[V r ,; 



-i \n, V r ] 



(179) 



Then, using the identity: [[A, B]_ , •] = [A, [B, »]_] + - 
[B, [A, •]+]_ the current superoperator anti-commutator 
can be decomposed as 



i (L n+ L v ' r - L V ' r L n+ ) 



(180) 



Here we introduced the cmiicommutator superoperator 
for the particle number (cf Eq. (f4"3")0 . 



L n+ = n,. 



(181) 



and decomposed the tunneling interactions into the junc- 
tion contributions, L v = ^2 r L v ' r . Importantly, the 
last term of Eq. (|180[) is irrelevant when inserted into 



Eq. fT77|) since TrTrL 



D R 

form of L Vr . We obtain 



due to the commutator 



(Hit) 



-t-TrL"+ [TrL v ' r p(t) 
2d \r pw 



(182) 



Integrating out of the reservoirs and collecting terms into 
irreducible blocks (cf. Sec. ([II B 3[l ) one now obtains 



(F)(z) = ±Tri/ l+ E r (z) 



D 



L(z) 



p(to) 



UTrL n+ Z r (z)p(z) 



(183) 



where E r is just that part of the irreducible self-energy 
E where the latest (leftmost) vertex is associated with 
reservoir r. We can thus decompose 



E(z)=5> r (*). 



(184) 



For the stationary current (F) = \\mt-t ->oo{I r )(t) = 
\\m z ^io — iz(I r )(z) we then obtain the central result of 
this subsection 



(F) = -iiTrL ,i+ E r (zO) /9 



D 



(185) 



where p is the stationary density operator (cf. 
Eq. CL52])). 

The first advantage of the Eq. f|l 85[) is that it allows 
one to explicitly see that the current is always zero at zero 
bias. For any number of electrodes, at zero voltage bias 
all electrochemical potentials are equal, p, r = 0, implying 
that all partial self-energies Y7 are proportional to the 
total self-energy is proportional to this value: 

E r cx E(i0) (186) 

L 'v T 

We note that this proof depends on the assumption 
of reservoir frequency independent spectral densities, 
r r (w) = r r , which we make throughout this paper. We 
may furthermore add L to E(z0) in Eq. (|185[) and express 
the current in the effective Liouvillian L(z) — L + E(z) 
since by local charge conservation, Eq. (f2"3")l , Tr L n+ Lu = 

Tr n[H, •] = Tr [n, H\» = when acting on any dot oper- 
ator. We see that there is no stationary current at zero 
applied bias, 



cx -iJrU 

D 



~L{iO)p = 



(187) 



since the stationary state p is the zero super-eigenvector 
of L(i0) (cf. Eq. (|37p). Importantly the above proof 
holds no matter what approximations one makes for the 
calculation of the self-energy E as long as all reservoirs 
are treated in the same way. This applies to perturbation 
theory up to any finite order, as well as to the RT-RG 
approach that we set up in Sec. (|IIIj) . 

A second advantage of Eq. (|185|) is that we can easily 
directly relate the current to only a few supermatrix ele- 
ments of the effective Liouvillian L(z) at zero frequency 
z = iO in the basis introduced in Eq. (|137p . The dot 
trace combined with the action of L n+ can be expressed 
in the dual supervectors of (|103[) and (|107[) 



1 



Tr(L 



2 D 



Tr(nu) 

D 



V2{T \+V2{Z L 



(188) 



Eq. (|85|) implies probability conservation (|69l) . but also 
more strongly that (Z^|E r (zO) 

this we obtain 



2TrE r (iO) = 0. From 

D 



(F) = -V2i(T \V r \p) 



(189) 



Clearly, the partial self-energy E r has the same general 
form as Eq. (| 137[) and we distinguish its parameters with 
an additional reservoir superscript r 



iE r = 

2T r \Z R )(Z R \ 



(190) 



C\Zr.)(Z l 



\Zr)$\ 



M r a \T a )(T a \+E> a \S a )(S a 



0){Zl\ 



\ ' F r 
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Inserting this form and the explicit expression for the 
stationary state Eq. (|152p into Eq. (|189[) we obtain the 
final explicit expression for the average stationary current 
in terms of the self-energy expansion coefficients: 

1 



(n = ^= (T \r)-(To\ec L \i>) 



(191) 



We emphasize that this equation is exact, given that the 
coefficients %j) r a and £ r of the partial self-energies E r are 
known, from which \jj a = ^ r i\f a and £ = J2 r £ r also 
follow. We can thus calculate the current very easily 
if we perform all calculations for fixed reservoir index r 
at the latest (leftmost) vertex and sum over r to obtain 
Eq. (|184|) . Finally, we note that using Eq. (|19ip . one can 
check explicitly that when imposing particle-hole symme- 
try on the expansion coefficients the current vanishes in 
agreement with the result obtained above (which makes 
explicit use of /il = expanding E r as Eq. (|173l) . 

one finds that \ip r ) and £ r have no components involving 

\To). 



III. REAL-TIME RENORMALIZATION GROUP 

In this section we set up the renormalization group 
calculation of the effective Liouvillian L(z). It is based 
on the real-time perturbative expansion (RT-RG) for the 
"finite-temperature" self-energy £ in the causal represen- 
tation introduced in the previous section, cf. Sec. pi B 41) . 
The procedure is to calculate L(z) by introducing an RG 
flow of the "infinite-temperature" Liouvillian L and the 
vertex G as function of a decreasing energy scale cutoff 
A, with the initial condition that L = E + L (thereby 
including the vertex of the type G) and G the vertex of 
perturbation theory. There are a number of motivations 
for performing such an RG treatment of the perturbative 
series for E. 

• First of all, treating E perturbatively in G, while in- 
finite orders of G have already been resummed into 
E, would amount to an inconsistent treatment. We 
note that even in the non-interacting case U = 
we already need to do a renormalizcd perturba- 
tion theory (|79|) up to 2-loop terms to recover the 
exact result all quantities (i.e., not just the cur- 
rent). For strong interaction U such corrections 
become of increasing importance since the broad- 
ening in E of the various excitations of the QD is 
uniformly ~ T i.e. it is energy independent. At 
finite and especially at low T the strong interac- 
tion U leads to non-trivial modifications of these 
widths due to quantum fluctuations, resulting in 
different lifetime broadening for single- (SET) and 
two-electron excitations (COT), with a non-trivial 
voltage dependence. Moreover, the bare perturba- 
tion theory breaks down at these resonances and 
T = and in the renormalized perturbation theory 
(fT9")l . the low energy cut-off T is replaced by the 
imaginary parts ~ T in L. 



• Secondly, the resonance due to the Kondo effect 
causes even the renormalizcd perturbation theory 
to break down and 3-loop corrections result in the 
enhancement of Kondo exchange processes. These 
have been studied using the RT-RG extensively but 
only after a Schricffer- Wolff transformation from 
the Anderson model (effective Kondo-model)^!. In 
the regimes of large applied bias voltage and mag- 
netic field such 3-loop corrections can be neglected 
due to the dephasing of the Kondo correlations. 
This is the regime of interest in this work and we 
expect that the 1- and 2-loop RG approach to give 
a good first approximation here that deserves to be 
studied. 

• In order to go beyond the study^i by studying 
Kondo effect in the Anderson model at least a 3- 
loop treatment is necessary. The renormalization 
of the 1- and 2-loop terms that we study here will 
then couple to the 3-loop terms and still play a im- 
portant role. Therefore, our study of the 2-loop 
RG provides an important starting point for such 
a much more involved study that in particular will 
address the low bias and low magnetic field regime. 



A. Flow of Keldysh contractions - continuous RG 

In general, RG approaches to transport aim to elimi- 
nate the effect of reservoirs states, starting at high ener- 
gies, by incorporating it into a redefinition of the system 
parameters. Typically one eliminates the states them- 
selves. Here, in contrast, we successively suppress the 
occupations of the states using an RG procedure while 
keeping the states. Before we specify the details and 
advantages of this cut-off scheme, we first outline the 
main idea of the functional renormalization group ap- 
proach when applied to the real-time perturbation se- 
ries (|79[) . By our causal reformulation of the perturba- 
tion theory (cf. Sec. pi B 2[) ) it is clear that all the in- 
formation about the occupations of the reservoir states 
is contained in the Keldysh components of its correla- 
tion functions, i.e., in the 7 contraction. Therefore we 
introduce a cut-off dependent contraction function 7a 
which monotonously flows from the initial, full contrac- 
tion function 7a|a=oo = 7 given by Eq. (|64[) to the trivial 
final function 7a|a=o — where all occupations are sup- 
pressed. During this flow we demand that the effective 
Liouvillian remains invariant: for every value of the cut- 
off parameter A, 



L(z)=Z + E({7,Z,G}) 

= Z a + E({7 A ,Z a ,Ga}) 
= La I 



IA=0 ' 



(192) 
(193) 
(194) 



Thus L(z) has the same functional dependence on the 
contractions 7a, the Liouvillian La, and vertices Ga- The 
latter two now acquire a A dependence to maintain in- 
variance. The same diagram rules are thus valid for any 
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value of A. As a result, at the end of the flow where 
7a|a=o = 0, the effective Liouvillian is given simply by 
La|a=o- The information about the reservoir degrees of 
freedom, previously incorporated into the self energy E, 
has now been incorporated fully into the dot Liouvillian, 
cf. Eq. (fl"94l) . 

The RG flow is generated by making an infinitesi- 
mal change dA of the cut-off, resulting in a infinitesimal 
change dj\ sa (dj\/dA)dA of the Kcldysh contraction 
function. In the perturbation series at scale A one splits 
up the contraction function as 7a = TA-dA+^TA and col- 
lects all terms in a perturbation series of the same form 
but containing only 7A— dA contractions. In this process 
the terms containing one infinitesimal contraction d"f\ 
can be identified with rcnormalizations dL of the Liou- 
villian, dG of the vertices and newly generated higher- 
order vertices with more than one leg. The process is 
illustrated in Fig. |4j The Liouvillian Za-<ja = L\ — cLLa, 
and vertices (?A-dA = G\ — dG\ of the new perturba- 
tion series are then all defined for the new, lower cut-off 
scale A — dA. We obtain differential equations for these 
quantities describing their renormalization as one con- 
tinuously reduces the cut-off A. These are the real-time 
renormalization group (RT-RG) equations. 

A key requirement in setting up this continuous RG is 
that for any A the zero-eigen vector of L A does not appear 
in the resolvents (z + X — Z A ) _1 , to avoid divergences 
as function of the frequencies. The RG thus has to be 
formulated such that the property Eq. ([55)) of the vertices 
in the causal representation is preserved. This can be 
shown to be the case, see Eq. (|214|) below. 

The final key point is the choice of a cut-off dependent 
distribution function in the contraction 7. The numerical 
integration of the RG equations is more stable when we 
introduce a contraction function with a cut-off on the 
imaginary frequency axis: 



(z-La-Sa)- 1 










+ 
+ 




FIG. 4: Renormalization group transformation with an 
infinitesimal change dA > of the flow parameter, A — > 
A — dA. The physical restriction is that the effective Liou- 
villian, L(z) — La + Ea, or equivalently, the density operator 
propagator (z — L(z))~ 1 remains unchanged. In lines 1-2 we 
start from the perturbation theory at scale A and split up 
the contraction function 7a (black curved lines) into the con- 
traction function with reduced flow parameter (7A-dA, blue 
curved line), which should appear in the renormalized pertur- 
bation series and the change the contraction function (~d"/A, 
red curve line). Next, in lines 3-5 the latter terms of linear 
order in d7A are collected into 1- and 2-loop renormalizations 
dZ/A of the Liouvillian in the resolvents and 1-loop renormal- 
ization dG?A of the vertices. Finally, the perturbation series 
is rewritten in terms of new Liouvillian Za-(JA and vertices 
Gh-dh defined on the new scale A — dA (indicated by blue). 
This transformation is exact if one also accounts for the gen- 
eration of higher order vertices^, which we, however, neglect 
here (they are not drawn). We do account for the renormaliza- 
tion of the original vertices, i.e., of single-charge fluctuations. 



7l2,A(W =«12-T> = —7 



through the function 



9 T (w) 



0(w), M>7rT 

U < 7TT 



ItxT ' 



(195) the contraction function (|195|) reduces to the simple form: 



(196) 



r r 

7i2,a(w = <5l2— / 



dJ (198) 

A rjuj — fi r — ibj 



See Ref. Fig for a detailed discussion^ 



where 0(w) is the theta- function and 
J = (21 + 1)txT 



(197) 



is the Z-th Matsubara frequency (I = 0, 1, 2, . . .). In the 
limit A — > +00 we recover the partial fraction expansion 
of the meromorphic function (r/27r) tanh(w/2T) as re- 
quired. Imposing this cut-off in Matsubara space leads to 
a suppression of the tails of 7 A (?7w) on the real frequency 
axis as A — > 0, rather than sharp truncation above fre- 
quency A. This implements the suppression of contribu- 
tions from states above energy scale A. From hereon we 
will consider the zero-temperature limit T — > for which 



B. RG in frequency space 

1. N on- equilibrium Matsubara representation and frequency 
dependence 

To formulate the RG equations we need a more com- 
pact notation for the various frequencies. Since the 
contraction functions 7 A depend on -quo = x — rjfi we 
re-express the reservoir energies in the resolvents in E, 
Eq. (|79)) , in explicit calculations as 



E-Xi = E- xx. 



E 



l...n 



(199) 
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(a) dL A /dA = 



E 



(b) dG A /dA 



dL\(E, lo) 



E\ 





6Eb 











dG\{E, uj, loi) 

FIG. 5: RG-equations and frequency dependence generated 
by the renormalization transformation of the diagrammatic 
perturbation theory in Fig. SI (using the same red and blue 
color), (a) Liouvillian renormalization by 1- and 2-loop cor- 
rections. The example diagram illustrates that the correction 
depends on the external frequency E of the diagram (Laplace 
variable) and on the sum of the reservoir frequencies ui run- 
ning over the dL/dA- block (marked yellow). These frequen- 
cies are read off at the vertical cut (green dashed line) to 
the left of this block, (b) Vertex renormalization by 1-loop 
corrections. The example diagram illustrates that the vertex 
correction depends additionally on frequency ui\ of the vertex 
leg- 



Here 1 , . . , n are the multiindiccs of the contractions going 
over diagram segment i. The frequencies are now taken 
relative to the electrochemical potentials and we write 
their sums as repeated multiindiccs: for k = r\k^kTkLOk^ 



Xk 


= Ul k + Pr k , 


Xl. 


.n = Xi + 




(200) 


L0k 








...+ U) n , 


(201) 


Mr 


= Vk^r k 


pi. 


.n = Ml + 


■■■fin- 


(202) 



Similarly, we express the dot energies relative to these 
chemical potentials as 



El 



^ 

i—l...n 



(203) 



A key advantage of the cut-off parametrization (|195|) 
is that it allows us to analytically perform all the in- 
tegrations over the reservoir frequencies uj in Eq. (|79j) 
by closing each integration contour in the complex lower 
half-plane and using the residual theorem. As a result, 
we can replace all integrations over real frequencies in the 
resolvents by summations over Matsubara frequencies ly- 
ing in the lower half-plane and Eq. (|199j) becomes 



El 



IL0\. 



(204) 



For finite T the sum of Uk runs over the positive discrete 
frequency values (2/^ + 1)tvT, Ik = 0, 1,2, .. which turns 
into an integral over positive u>k for T = 0. As always, 
we do not explicitly indicate this integral^ 
2. RG equations 

During the RG flow the Liouvillian develops a non- 
trivial dependence on both the real energy E of the QD 
and on the sum of the imaginary frequencies of the reser- 
voirs, iui, of all contractions that pass over the propaga- 
tor in a diagram. This is illustrated in Fig. EJa). We 
separately keep track of these frequency dependencies by 
writing 



L A (E,u) := L A (E + iuj), 



E + - L A (E,u) 



(205) 
(206) 



The renormalization of the vertex G\ introduces similar 
dependencies and an additional dependence on the reser- 
voir frequency u>\ of the vertex "leg", i.e., the contraction 
connecting it to another vertex, see Fig. [5jb): 



Gi iA (£',a;,a;i) := G x ,k{E + iw,iu{) 



(207) 



The formally exact, infinite hierarchy RT-RG equa- 
tions were derived in general form in Ref. I2H . Here we 
restrict ourselves to the 1- and 2-loop order approxima- 
tion for the Liouvillian and the limit of T — > 0: 



dL{ -f: ^ = i-G x {E, u, A)II(£i , lo + A)Gi (E u lj + A, - A) + L x (208) 
aA it ir^ 

Gi (E, u, A)n(£i , u + A)G 2 {E x , u + A, uj 2 )II{E 12 ,lo + A + lj 2 )G 2 (E 12 ,lu + A + cj 2 , -cj 2 )II(£?i , u + A)Gi(E u u + A, - A) 
This approximation requires that one also accounts for the renormalization of vertices to 1-loop order: 

ddiE uj,^) = _JiQ a ( E MjLtfy + a)G 1 (E 2 ,lo + A 1 lo 1 )U(E 12 ,lu + lo 1 + A)G 2 (E 12 ,lo + lu 1 + A, -A). (209) 
aA ir 



I 

Both here and below we leave the A dependence of the renormalized L and G and the summation / integration 
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over all internal indices implicit 

By construction the self-energy is obtained by 

directly integrating Eq. (|208j) : 

E(E + iuj) = L(E + iu)) - L = \ dA—^(E,Lu) (210) 

For the transport current we need reservoir- resolved parts 
S r (z) of this self-energy (cf. Eq. ([184])), which is obtained 
in the same way: 



f] r (E + ioj) = L r (E + iuj)-L r = I rfA-^- 



(211) 



where L r A \ A=oc = L r := XT_ is given by Eq. (J7T}. The 
RG-equations determining L\ and its associated vertex 
G r have the same structure as the RG equations from 
Eq. (|208j) - (|209j) and are simply obtained from the latter— 
by providing the left most vertex on the right hand sides 
of Eq. flU} and flU]) and the left hand side of Eq. ([305]) 
by a superscript r. We emphasize that ^ r L r — X = 
L-L and J2 r L A \ A = = S(z) = L(z) - 

The frequency dependence in the RG-cquations 
Eq. (|208p - (|209p in the 1 and 2-loop approximation needs 
to be carefully discussed. Before we turn to this in 
Sec. (|III C|) - (|IIID[) . we discuss important general prop- 
erties of the RG-equations in Sec. (|III B 3[) and their im- 
plications for the fermionic eigenvalues in Sec. QUI B 41) . 
Using these results, we can decouple some of the RG 
equations (see Sec. Cp ~) and show that in the 1-loop 
approximation we already obtain the exact solution for 
the current in the limit {7 = 0, even though in this limit 
2-loop corrections are non-zero (see Sec. (|III C 3|) ). 



Eq. (|214[) - (|215l) follow by assuming that they hold for a 
given A, and by (Z L \^G = ^G\Z R ) = , obtained 
acting with Eq. (|209|) on these vectors. Since Eq. (|214|l - 
(|215l) hold initially for A = oo the result follows for any 
A, E, u>, and uj\. 

Similarly, we now show that for any A, E, to, and u)\\ 



d 

dA 



L{E,w) 



n, 



0, 



(a+\SrL(E,u>) = 0, 



(216) 
(217) 



In our 2-loop approximation for dL/dA this follows from 
the property of the renormalizcd 1-leg vertices 



G(E, u},Wi)\a na .) oc \Z R ) or 
(a+ r |<5(£,w > wi) cx (Z L \ or 0. 



(218) 
(219) 



The proof of Eq. (fT2"g|) - (fT2"2"Jl can be extended to the 
renormalized vertices as follows. We start by observ- 
ing that the right hand sides of the RG equations (|208|l - 
([HJ| have the same structure as T,(z) (cf. Eq. ([75]) ). 
Assuming that Eq. (|218[) holds for a given scale A, the 
RG equation (|209[) implies that it is preserved under the 

flow: (a+jJ^-G = ■a A -G\ a v< r ) = °- Here we uscd tnat 
by Eq. (pT2) -([2l3 ]) both \Z R )_a,nd (Z L \ arc eigen super- 
vectors of the renormalized L for all h.,E,w,uj\. Since 
Eq. ([2T5 ]) -([2Tg ]l hold initially for A = oo, this then im- 
plies it holds for all A. From this Eq. ([216]) - ([217]) fol- 
low directly. The above proofs are readily extended to 
the infinite hierarchy of exact RT-RG equations for ver- 
tices with multiple legs, confirming that Eq. (|212l) . (12131) . 
(|217| and (|216|) hold exactly (and not just in our 2-loop 
approximation) as required by Eq. (|210[) . 



Exact eigenvectors 



Since according to Eq. (|210p the self-energy 
is obtained by integration over A, the RG equations 
(|208[) - (|209[) must preserve the exact properties derived 
in Sec. pi G 5 a[) from the causal structure of the pertur- 
bation theory, Eq. (fT9"|) . Indeed, an exact left and right 
eigen supervector of the renormalized L is given by 



(Z L \L(E,u) = 0, 
L(E,oj)\Z r ) = -i4T\Z R ), 



(212) 
(213) 



respectively. Since we have the initial condition at A = cxj 
(Z L \L = and L\Z R ) = -i4T\Z R ), we only need to 
show that (Zl\j^L = jj^L\Z R ) = respectively, to prove 
Eq. (|2"T2"|) - (|2"T5|) . We note that since (PTU]) is an exact 
relation, Eq. (|212[) - (|213[) must hold for the exact, infinite 
hierarchy of RT-RG equations (i.e., including all higher 
vertices generated during the RG-flow that we neglect 
here). In our 2-loop approximation this relation directly 
follows by letting Eq. (|208[) act on these vectors using 



(Z L \G(E,u),u) 1 ) = 
G(E,u,u 1 )\Z R ) = 0. 



(214) 
(215) 



4- Fermionic excitations 

a. Fermionic eigenvalues We can now pick up the 
discussion of Sec. pi Dp . Since the supermatrix structure 
of dL\/dA in the fermionic sector is preserved under 
the RG flow and is the same as that of we can 

now directly relate the coefficient AF~ a + introduced in 
Sec. (jITTT3l) using Eq. (I2TU1) : 



&F~+(E- 



dA- 



dF. 



■qa,A 



dA 



(E + iuj) (220) 



This coefficient determines the fermionic excitations at 
arbitrary complex frequency as given by Eq. (|164[) for the 
exact L{z) when the infinite hierarchy of RG-equations 
is used to compute the right hand side. We see that the 
A-dependent coefficient F~J~ A of La interpolates between 
the infinite temperature limit, where AF~ a + = 0, and the 
exact value Ai^ + of S(z) through Eq. I|220p as was an- 
ticipated in Sec. pi Dp . All renormalization effects enter 
into the fermionic excitations through the renormaliza- 
tion of the four complex coefficient F~ a + of L\. During 
this flow, the qualitative features of these excitations, dis- 
cussed in Fig. [3] may change. During the continuous RG, 



2G 



the complex parameters AF~^ + will grow from zero and 
modify both real and imaginary parts in Eq. (|164p . This 
happens only for the interacting Anderson model, U 7^ 0, 
since U multiplies these coefficients in Eq. (|164[) . This 
flow may include bifurcations as function of the flow pa- 
rameter A but for large enough U 3> 2r the excitations 
energies (real parts) remain non-degenerate. However, 
the general result Eq. (|169|) shows that during this non- 
trivial flow the average of the complex eigenvalues stays 
fixed for all frequencies. We conclude generally that the 
fermionic excitation energies and decay rates are renor- 
malized symmetrically with respect to the average values 
e + y + cr -j and 2T, respectively, for any complex fre- 
quency E + iuj with lo > 0. 

Finally, we note that the stability constraint discussed 
in Sec. piD 31) imposes a constraint on the RG flow: since 
at any stage of the RG flow the effective Liouvillian L{z) 
(|193| can be calculated from the pcrturbative expansion 
(|79jl . the imaginary parts of all non-zero eigenvalues of 
the renormalized L\ must be negative to avoid unphys- 
ical divergence of the time-dependent density operator. 
Such behavior would not go unnoticed in the RG, as it 
would lead to an instability in the flow since zero denom- 
inators would appear in the resolvents in Eq. (|229[) . This 
provides a simple criterion for the stability of the RG 
flow for the Anderson model that can be checked eas- 
ily in numerical approximations. Although in previous 
applications of the RT-RG no instabilities have been re- 
ported and in the present study none were encountered 
either, the general conditions for stability are currently 
not known. 

b. Fermionic supermatrix elements The properties 
Eq. (|218j) - (|219|) strongly restrict the fermionic matrix el- 
ements of resolvents II that can appear in the RG equa- 
tions. We will see that implies that quite generally the 
RG equations decouple into smaller sets of equations (see 
Sec. ()III C[) ) and that important simplifications arise in 
the U — limit (see Sec. ()III C 3jl ). These simplifications 
arise since in general on the right hand side of RG equa- 
tions such as Eq. (|208[) and (1209j) the resolvent II always 
appear sandwiched between pairs of G vertices (all are 
renormalized quantities but A is not written here). 

• In matrix elements of terms with only one resol- 
vent, (A|GilG|F) there are no restrictions only 
if X = Z R ,Y = Zl,. Indeed, upon inserting 
the completeness relation 1 = J2 V a v I^o-X^ctI + 
(bosonic terms) left and right of Id, we see 
that according to Eq. (|218[) - (|219[) all intermediate 
fermionic supermatrix elements contribute. How- 
ever, when the basis supervectors X, Y involve one 
zero eigen-supervector only certain matrix elements 
contribute: for v = ± 

(a; a \n a \a- a ), X = Z R ,Y^Z L (221) 
(a+ a \U a \a; a ), X ? Z R ,Y = Z L , , (222) 

whereas, if there are no zero eigen-supervector 



among X,Y only one factor is possible: 
(a+ \n a \a~ a ), X^Z R ,Y^Z L . (223) 

• In terms with n > 3 resolvents 
(X\GUG ...GUG\Y) Eq. ([22H and (f223|) ap- 
ply to the leftmost "boundary" resolvent. 
Otherwise the expression vanishes for any A 
since- by Eq. (X\GU\a + )(a + \GHG . . . ex 

(X\GU\a+){Z L \UG . . . = 0_where we used that 
(Zl\ an exact eigenvector of L and thereby of n by 
Eq. (|212|) and a zero eigenvector of G by Eq. ()214|) . 



• Similarly, in such terms Eq. (f222]l and ([223]) ap- 
ply to the rightmost "boundary" resolvents since 
...GnGn\a-)(a-\G\Y) = by Eq. ([2T5|l . (^131) . 
and (12151) . 

• Finally, in such terms with n > 3 resolvents, 
the resolvents which are not at the boundary 
can only contribute with fermionic matrix element 
(a+jn^la^), irrespective of X and Y: this factor 
must always occur at least n — 2 > 1 times. 



C. 1 loop RG equations 

1. Frequency dependence 

Since our goal is to calculate the stationary state 
from the effective Liouvillian L(z) = L(E,ui)\ A _ at 
z = E + iuj = iO, we first consider the RG equation 
for this quantity in the 1-loop approximation and at fre- 
quency E = 



dL°(0) 
dA 



z-G°n°(^,A)G°. 

7T 



(224) 



Here the superscript indicates that we also evalu- 
ate the Liouvillian at zero reservoir frequency, ui = 0: 
L°(E) := L(E, 0) Similarly, in Eq. ([2lM]) we approximate 
the vertices by their initial values, as given by Eq. (|1 19[) 
and (fT2T|) or Eq. (j5Tj). 



Gi(0,0,A)«Gi(/zi,A,-A)«G? 



(225) 



neglecting their dependence on the dot frequency (Ei = 
/2i), the reservoir frequency (uj = A) and the vertex- 
leg frequency (uji = ±A). Such frequency dependencies 
arises only when accounting for the renormalization of 
the vertices: for small frequencies we can approximate in 
Eq. dp, 



dG 
dA 



r 

X 1 



G'' 



(226) 



on the right hand side G ~ G° ~ 1, giving G = 
G° + 0(r/A). In 1-loop order for L one must therefore 
consistently neglect the renormalization of G, Eq. (12251) . 
with respect to the log-corrections to the Liouvillian that 
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arise from Eq. f|224|) . This will be checked later on. The 
resolvent is likewise evaluated at uj = A. Note that in- 
dependence of the resolvent in first approximation, 



i 



(227) 



does not originate from the Liouvillian. 

In contrast, Eq. (|224p depends on QD frequency E 
in an important way: due to the finite bias voltage the 
renormalization of the zero i?-frequency Liouvillian cou- 
ples to the finite frequency Liouvillian E — > E\ = jli = 
appearing in the (|206[) on the right hand side. 
We therefore need to consider instead the RG equations 
at finite QD frequencies, 

= i~G\U\E u A)G? (228) 

CIA TT 

which lie on a discrete grid: E = kLp.L + kRfiR = 
rj(k L - k R )V/2 where k r = 0,1,2,..., r = L, R. In 
the numerical calculations we keep as many equations 
as required to make the solution converge with respect 
to the k r . This coupling of the RG flow of the Liou- 
villian at energies differing by multiples of the voltage 
arises because the Matsubara frequencies of the different 
reservoirs are shifted by different electrochemical poten- 
tials: this is typical feature of renormalization in a non- 
equilibrium system^. We discuss the effect of neglecting 
the QD energy ^-dependence in the RG equations in de- 
tail when we analyze the numerical results in Sec. (|IVBI) . 



we obtain, abbreviating 0f = E\ — \ k (E\), 



^(k 3 \L°(E)\ Ko ) = -(k 3 \M\k ), (230) 



where \ni) are elements of basis (|91I) - ((M|) and the super 
matrix elements are sums (ki.2 sums implicit) of factored 
contributions: 



(k 3 \M\ko) = ^(«3|G°| K2 )( Kl |G?| Ko )(^|ni| Kl ) 



The product of G° matrix elements gives a simple numer- 
ical factor or ±1 (see expansions (|119[) - (|121[) ). whereas 
the super matrix elements 



»(K2|nil«i) 



i 



A - iO\ 



(k 2 |^|«i) 



arise from the spectral decomposition of the resolvent, 
IT (Ei , A) = III ■ To explicitly sum over r\ (contained 
in the multiindcx 1 = 77, a, r) which enters the resolvents 
only through E\ = E — r]fi r := E vr we abbreviate for 
now A*(£i) := A, )r , P^E-l) := P vr and 



2. Explicit 1-loop RG-equations for the Liouvillian 

Inserting the spectral decomposition for L(E), c.f. 
Eq. (H121), 



dL°(E) _ T 1 



dA 



tt A-iQ* 



(229) 



" r A - iE w + i\\ )r 



(231) 



Expanding Eq. (|2"2"5|) in basis Eq. ([9"T]). Eq. we ob- 
tain: 



dL°(E) _ iT 
dA TT 



- (aI ff |lC;>^) + teJn£' 



Z R )(Z L \ + (al^llCr l«=a) " (a+tKT 



\Zr)( X o 



(aU\nl7\a^) (aX a \nl+°\aX„))\x»)(ZL\ (Ja+JII^KJ + (<*U\n°Z'\<*Z ff ))\x*)(x>\ 
(a+ a \Kr\a+ a ) + (at a \U a _-'\aZ a ))\S a )(S a \ - (aU^-V \a^)\T + )(T + \ (at tr \I&'\eC <r )\T--)(T-\ 



{X^\^\x^ + (T4n^\T-)-(S a \T^\S a ))\aZ a )(a± a \ 



(232) 



I 

Here we leave implicit the summation over a and r, as and oi„ a subspaces (cf. Eq. (|162[1 ). The first two terms 
well as the summation over the two eigenvalues in the x m the equation don not contribute to the calculation the 
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remaining terms of the effective Liouvillian or the trans- 
port current but are written here for completeness. Due 
to (Zl\G = (cf. discussion of Eq. (|85[l) the zero eigen 
projector P Zl does not appear in Eq. (|232p . ensuring 
that none of the resolvents can diverge during the RG- 
flow. 

Using the relations K(A - i(E - L(E)))~ 1 K = (A - 
i(-E* - L(-E*)))-\ cf. Eq. ([70]), (HT7J and CUE]), and 



expanding resolvents into eigen-projectors Eq. (|231l) we 
obtain the explicit RG equations for the Liouvillian ex- 
pansion coefficients. We first have a set of equations for 
10 complex coefficients on an infinite, discrete grid of fre- 
quencies E 



d£a,a(E) = .T( (a + _ s \P"- r s \aZ s ) ^ H*\ p +? 
dA \ I A - i(E_ r - \ a _~ s 



dE a (E) 
dA 

dM+(E) 
dA 

dM- (E) 
dA 

dFt~{E) 
dA 

dF+~(E) 
dA 



t / {a\ a \ptrw- +a ) 

\ \ A - i(E-\ 



= —i 



.r (aUP- 



A-i{E +r -\ a + X'\ 
A - i{E +r - X 



7T A - i (E- r - A 
r «JP^>-J _ dMl(-E*) 



—i- 



—i 



i{E +r -X a + --) 
-i(E- r -\* r ) ' a -*(£?. 



dA 

(t + \ p+;\t + ) 

'+r - A + + ) 
(T_|P^|T_) 



7T I A - i(E +r 



dEU-E*) 



dA 



(S a \P*°\S a ) 

A - i(E- r - X-r) j 



A - i(E_ r - \ T r) A-i(E^ 




d(F+r(-E*)Y 
dA 



(233) 
(234) 
(235) 
(236) 
(237) 
(238) 



Importantly, the eigen projectors of the Liouvillian L , 
PI = P l (E — ri(i r ) with eigenvalues X^ = X l (E — i]fi r ) (cf. 
Eq. (|149[0 depend on the frequency E~i]ii r . The explicit 
expressions for the projector matrix elements on the right 
hand side are given in Eq. (|151|) . (|153|) . (|154p(|156j) and 
(|160[) and involve only the 10 coefficients appearing on 
the left hand side. Eq. (|233|) - (|238[) thus form a closed set 
of equations. This derives from the fact that the eigen- 
projectors of La that involve the zero eigen supervectors 
Zl and Zr drop out on the right hand side by Eq. (|162l) . 
Note that the fermionic matrix elements in the equations 
for the coefficients of bosonic terms that do not involve 
a Zl or Zr (£ CT(T ', E a and M v ) display the restriction 

msa . 



right hand side of Eq. (|2l?3"|) - (j2"3"8"|) . Their RG-equations 
dME) 



dA 



(239) 



F / { a ± 9 \P^'\a± s ) {a\ a \PlrH 



7T I A — i(E- r - X_~*) A - i(E +r - A++ CT ) 



dC(E) 
dA 



V / (aZ s \P«r\aI s ) 
\ { A-i(E- r - X a _7 ) 



dME) 
dA 



(240) 

A- l (E +r -X a + r) / 
(241) 



r / (aZjPl 



A - i(E^ r - X°L~° ) A - i(E +r - Xll, 9 ) 



The following 5 complex coefficients do not appear in 
the eigenvalues and projector matrix elements on the 



are therefore not required for the solution of Eq. (|233|) - 
(|238[) . but do depend on the solution of the latter. How- 
ever, for the calculation of the current only ip a is re- 
quired. In contrast, the coefficients Q and 4> a and are not 
required. However, these coefficients do renormalize and 
are required if one wishes to calculate, e.g., the stationary 
density matrix (|152[) . 

The remaining coefficients do not flow under the 1- 
loop RG, and remain at their initial values. For the two 
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diagonal matrix elements of the coefficient matrix £ in 
the bosonic sector we have 



d£,a 



dA 



<E) = 0, 



(242) 



which is valid only within the present f-loop approxima- 
tion. In contrast, for the remaining 16 fermionic coeffi- 
cients we have in general (e.g., also in 2-loop order) 



dF +± 



dF ±_ 



0, 



(243) 



due to the causal structure (cf. Sec. (jll D 3|) and 
Sec. (|IIIB4p ). We furthermore note that Eq. (|233|) - ((24Tj) 
explicitly satisfy the Hermicity conditions Eq. (|145[) and 
(|148[) . Moreover, the proper transformation under charge 
and spin-rotations is explicitly guaranteed by our use of 
irreducible tensor supcropcrators, cf. Sec. (jll C 5 b|) . 

Finally, for calculation of the current at a specific elec- 
trode r = ± we need RG equations for the coefficients of 
its self-energy component L r , cf. Eq. (|21ip . These are 
simply obtained from the above equations by (i) giving 
all coefficients a superscript r and (ii) suppressing the 
summation over r contained in the multiindcx 1 on the 
right hand side, i.e. setting 1 = r], a, +,u>. 

Before we proceed to calculate the 2-loop corrections 
to Eq. (|224[) we first show that already in the above 1- 
loop approximation we obtain the exact solution for the 
current in the limit U = 0. This is important since it 
demonstrates that for the current both the 2-loop cor- 
rections to L, the 1-loop corrections to the vertex G, and 
the uj frequency dependence that we neglected here, are 
intimately connected with interaction effects. We note, 
however, that for U = there are non-zero 2-loop cor- 
rections which do not affect the current^. 



3. Non-interacting case U = 0; exact solution 

Without local interaction, U = 0, the Hamiltonian 
Eq. (Ilip is quadratic in the fermionic operators and the 
non-equilibrium Anderson model can be solved exactly 
in this limit. A solution using the real-time approach 
was reported in Ref. I5l"ll52l We now show (i) how us- 
ing superfermionic algebra one can obtain this solution 
within RT-RG framework and (ii) that within the 1-loop, 
frequency-independent approximation (|224p this result is 
recovered upon careful inspection. 

a. Exact current In general, to calculate the cur- 
rent according to Eq. Q19ip . we need the elements £ and 
i\) a of the bosonic part of the effective Liouvillian L\. 
From these we can then easily find the other required 
coefficients £ r and ip^ of L r A which we do at the end. 
The coefficients £ and ip a in turn require the solution of 
Eq. ([2"3"3>([2"3"g|) which we first discuss. We show then 
that higher order correction as well as frequency correc- 
tions which we neglected in deriving Eq. (|233[) - (I238|) have 
no influence on the stationary current. We start by not- 
ing that on the right hand side of RG Eq. ([2l?3> (|2ll6l for 



the bosonic sector only resolvent matrix elements appear 
with intermediate fermionic states. We therefore first 
calculate the eigenvalues of the fermionic projectors to 
which these coefficients couple. From Eq. (|164p it follows 
that for U = 

\ a *<"± = rje + - -2iT±iT, (244) 
P"*"* = \ (< ± ^< T V^f^^-^j (245) 

Importantly: (i) Since U = the eigenvalues are inde- 
pendent of AF~ a + , i.e. they are not renormalized (cf. 
Sec. (|IID 3P ) and therefore do not acquired a frequency 
dependence (ii) The right hand side of the RG equation of 
any bosonic supcropcrator that is relevant for the current 
(i.e., excluding C,,<p), with the exception ijj a by the gen- 
eral property (|222p contain the off-diagonal super matrix 
elements 



(a+ |IC>- CT ) = 



(246) 



as a factor. It vanishes for U = by Eq. (|245|) . This 
implies in fact that that these coefficients do not renor- 
malize in any higher loop order since such terms on the 
right hand side of Eq. (|208p always contain this factor at 
least once by Eq. f|223|) . (iii) The renormalization of the 
quantities ip a at E = involve fermionic virtual states 
with the simple factors 



r 



i 

al A - irV/2 + i\l + r ° (rV/2) l ~ 
1 



A + T + i{e„ - rV/2) 



14.) (247) 
(248) 



instead, where e a :— e + aB/2. Importantly, for U = 
these matrix elements are independent of Af^ + . Higher 
loop corrections for ip a vanish since they contain the fac- 
tor (pM6"P at least once by Eq. f2"2"Tj) . Therefore Eq. ([2"3"9"P 
is the exact RG equation for ip a at E = and for U = 0: 



#q(0) 

dA 



E 



2r, 

2T 



£-M«tr|m£-|a+,) 

e CT - rV/2 



7T (A + T) 2 + (e CT - rV/2f 



(249) 
(250) 



where we again used the A'-conjugation properties 
Eq. jTDJ), pTTP and ([TTgp . With the initial value 
Vv,A=oo = (cf. Eq. (|135p .'l we obtain 



Vv(0) = — — arctan ' 



Leaving out the summation over the electrode r in the 
above calculation, wc obtain the coefficients of the self- 
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energy component L^|a=o = ^'"(0), required for the cur- 
rent (cf Eq. CDS]) ), 



= iE(#-#)' 



(251) 



using 775 (T | = 5X)(x<r| and £ r = £/2 giving finally 

(X 



(-0 = r ~ 



r.cr— ± 



r 
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r 



and for the non-linear differential conductance: 

r 2 



di_ 

dV Air 



- y 



r,(T— ± 



r 2 + (e a + rV/2Y 



(252) 



(253) 



in our units e = 1, ?i = 1 Restoring CGS units the current 
and conductance prcfactors become — > T§ — > 
giving in linear response a conductance of e 2 /h per spin 
channel. The importance of recovering this exact result 
for U — is that already at this level of approximation 
our the RT-RG approach captures correctly the weak in- 
teraction limit U <C r. while treating the tunneling non- 
perturbatively in L. Moreover, it shows that the 2-loop 
corrections to the Liouvillian and 1-loop corrections to 
the vertices which affect the stationary current are gen- 
erated by the Coulomb interaction. In general, however, 
the non- interacting limit requires a 2-loop treatment. 



D. 2 loop RG equations 



Vertex frequency dependence 



As concluded in Sec. pil C ip in the 2-loop approxi- 
mation for the Liouvillian we should consider the vertex 
renormalization and w-dependence of both G and L. In- 
deed, we find below that these effects are comparable 
and involve important cancellations. By systematically 
expanding about the frequency-independent bare vertex 
G , we can incorporate the vertex corrections into a sin- 
gle effective equation for the 2-loop Liouvillian, Eq. (|272[) 
below. We proceed in three steps: 



Step 1 Starting point is the 1-loop approximation de- 
fined by Eq. (|228p for any E. This we use to calculate 
a first approximation for the w-dependence of both the 
Liouvillian and the propagator. We expand 

L(E,w) « L°(E) + L 1 (E,lo) (254) 
IL(E,u) « U°(E,uj) + n 1 ( J B,w) (255) 

The 1-loop equation accounting for the leading frequency 
dependence is obtained by setting G sa G° and H(E, u>) sa 
U°(E,u)) the 1-loop part of Eq. (gOg) : 

dL(E u) = Jgjojp^ A + (256) 
dA 7T 
Subtracting Eq. (|228|) . we obtain 



dLl f k u) « £g\ (n°(Bi, « + A) - n°(£ , i, A)) G\. 



Shifting the integration variable in the II(i?i, A) term we 
obtain (in wide band limit assumed throughout): 



L\E,w) 



i—GiU°(E 1 ,u} + A)Gj (257) 

7T 



A-(, 



Eq. (|257|) docs not need to be calculated further since it 
cancels out below. Note that the correction vanishes at 
zero frequency, L 1 (S,0) = 0, for all E as required. Ex- 
panding the full resolvents Eq. (|206p with the approxima- 
tion Eq. (|254[) to the first order in Z 1 (i?, w), we obtain: 

n 1 ^,^) = n°(£>)I 1 (£>)II (.E,a;) (258) 

Step 2 In a similar way, we now calculate the ui cor- 
rections in the leading 1-loop order for the vertices: 



Gi(-E,w,wi) « G° + G\(E,uj,uj 1 ) 



(259) 



Keeping on the right hand side of Eq. (|209l) only the 
leading term we obtain with the same approximations as 
above 



Gl(E,u,U!) = -i- I dA'G° 2 U°(E 2 , 

D 



A')G° 1 n°(E 12 ,uj + Lj 1 + A')G° 2 - 



(260) 



We stress that the argument under the integral depends through the cut-off dependence of La' (G° is the bare 
on A' both through the explicit arguments as well as 
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vertex). Restoring the latter explicitly, 



cj + A' — iL° A , (E) 



(261) 



Since this is the only point where this is important we 
stick with the implicit notation. 

Step 3 Using the expansions for the resolvents 
Eq. (|258p and vertices Eq. (|260p we can now calculate 
an approximation to the right hand side of Eq. (|208p . 
keeping the leading order frequency corrections: 



dL 
dA 



di (lloop) d£(2loop) di(^rtex) 



dA dA dA 

all written at frequencies E and w, where 

dU 1 lo °p) . r 



(262) 



i— 

n 



G 1 U (E u A)G 1 + G° 1 Il 1 (E 1 ,A)G 1 

(263) 



are the terms appearing from expansion of 1-loop dia- 
gram in II 1 , 

J j (2 loop) r 2 

^— = L -G° 1 U (E 1 ,A)G° 2 (264) 
x Il (E 12 , A + Lj 2 )GllL (E 1 ,A)G° 1 
is 2-loop term with bare vertices, and 

A J (vertex) p 

= i-G\{E,Q,A)V°{E x ,A)G\ (265) 
dA it 

+i-G° 1 U°(E 1 ,A)GUE u A,-A) (266) 

7T 

are the terms appearing from the expansion of G in G 1 in 
the 1-loop diagram. Inserting Eq. (|257|) into the 1-loop 
Liouvillian frequency correction (last term in Eq. (|263p ) , 
we see that it exactly cancels the 2-loop term Eq. (|264l) . 
Neither term therefore needs to be calculated, simplifying 
the approach to a great extent. 

The above holds for any E and u>: integrating the 
Eq. (|262p at u> = using the above calculated right hand 
side we obtain a new approximation for L (E), improv- 
ing over our initial approximation based on the 1-loop 
equation Eq. (|228|) . In principle, with this steps 1 and 2 
should be repeated, resulting in corrections of higher or- 
ders which we neglect. We thus equate dL/dA ps Lr /dA 
on the left hand side of Eq. (I262p . We obtain a central 
result of the this section: an effective 2-loop RG equation 
for the Liouvillian at uj = 



dL°(E) r 



dA 



~G?n°(£i,A)G? 
i-Gi(JS,0,A)n°(Bi,A)G? 

r 



7T 



G^TSP(Ei,A)G\(Ei,A > -A) (267) 



Notably, due to the cancellation the entire leading reser- 
voir (w) frequency dependence comes from the vertex cor- 
rections Eq. (|260p . This single equation yields a signifi- 
cant simplification over the coupled equations Eq. (|208p 



and Eq. (|209|) . As expected, the elimination of Eq. (|209|) 
makes Eq. (|267p an intcgro-diffcrcntial equation for L. 
However, it can be converted into a differential RG equa- 
tion by analytically performing the integral in Eq. (|260p , 
see Sec. ([Ill D 2p . Furthermore, since we can work with 
the bare vertex super operators and we can make use 
of their simple anti-commutation relations Eq. (|53[) , that 
arc not preserved under the RG (in contrast to other 
useful properties of the vertex, see App. [Gj. The QD fre- 
quency (E) dependence in Eq. (|267p is of the same type as 
for the 1-loop equations Eq. (|232p : the RG equation for 
L°(0) depends on L°(fii) = L°(r) 1 r 1 V/2), etc. It there- 
fore has to be solved in the same way by including multi- 
ple Matsubara axes and converging the energy-hierarchy 
of equations, see Sec. (jlll C l|l and Sec. (|IV|) . Finally, we 
also note that one can indeed neglect the frequency de- 
pendence generated by vertex renormalization since it is 
indeed small, as we assumed in our derivation of the 1- 
loop equations. This can be seen if one substitutes the 
calculated correction Eq. (|260p into Eq. (|209p . Here we 
anticipate the projector expansion of G 1 , Eq. (|27ip : it 
is seen that G 1 is a well behaved function of the cut-off 
and frequency, decaying at small A, and it generates only 
small corrections in agreement with our approximation 
Eq. AMU) . 



2. Explicit 2-loop RG equations for the Liouvillian 



To obtain a differential RG equation from Eq. (|267l) . 
the integration Eq. (|260p needs to be performed. This 
complicated by the implicit A' dependence in the prop- 
agators pointed out with Eq. (|26ip . We now make an 
adiabatic approximation by expanding only this depen- 
dence about A' = A, i.e. we substitute 



La 1 ~ L\ 



(268) 



in Eq. (|26ip and neglect corrections ~ dL A i / dA which, 
by the RG equation Eq. (|267p . are of higher order and 
should therefore be neglected. To preserve the compact 
form of the equations we define 



PL n = P i (E 1 ... n ) 



(269) 
(270) 



where P 1 (E) and \ % {E) are eigenprojectors and eigen- 
values from Eq. (|149[) at cut-off A (not A') and we will 
implicitly sum over all appearing eigenvalue labels k 
below. To perform the integral we insert the projector ex- 
pansion Eq. (|149p of L evaluated at A under the integral 
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and obtain the explicit u;-depcndcnt vertex correction: 
G\{E,w,Wi) = (271) 



G° 2 P\E 2 )G\Pi{E 12 )Gl 

n I (a' + uj- «e|) (a' + oj + ui - ieQ 



dA' 



,T G° 2 P'(E 2 )G° 1 P^ {E 1 2 )Gl 

n LJt-i (e{ 2 - e 



In 



A + uj-i<d l 2 
A + uj + uji — iQ\ 2 



Combining the rest of Eq. (|263]) with Eq. ([2651) we obtain: 

dL°{E) r i 



— rG?P fc G? 

;ofc ill 



dA 

r 2 



7r A — iQf 
1 1 



(272) 



7r A — iOf A 



12 



■In 



2A - iQ{, 
A - iOi, 



^1^1 ( - r 2^12 Lr 1^2 ( - r 2 + ( - r 2 i 2 ( - r lM2 ( - r 2^1 °I 



pi stO pi ,^0 pA; ^<0 



This explicit result is a central result of this paper. 

The explicit evaluation of Eq. (|272[) for the Anderson 
model is required for our numerical implementation, but 
also allows us to draw some general conclusions about 
the 2-loop (and higher) corrections. The 1-loop part is 
that given by Eq. (|228[) and we proceed analogously for 
the 2-loop part: 



dL°(E) 



dA 



2 loop 



-i—( K7 \M\K )\K 7 )( Ko \ (273) 



where the super matrix elements are factored as follows: 
(k 7 \M\ Ko )= (274) 

]T («7|G;|K6)(KB|G§|«4)(«3|G?|K2)(Kl|G§|«0) 



K 6 ..Kl 



x Af ( {k 6 \P?\k 5 ), (k 4 \P( 2 \k 3 ), ( K2 |P 2 >i) 



The product of G° matrix elements gives a simple nu- 
merical factor, whereas the remaining part, 

7V((«6|i^|K5),(«4|i?al«3),(K2|i^|«l)) - (275) 

5(0j , @l Q{ 2 ) (k 6 |P* I « 6 ) (/s 4 | f? 2 1 «a) («a 

contains the product of the non-trivial projector matrix 
elements and the propagator factors 



I A-i0i 



(A - i9f) (A - * (6{ 2 - 9^ 

(ef o e|) (276) 



The argument of this function is constructed by formally 
putting the variables containing the eigenvalues of the 



three projectors in the argument M into the correspond- 
ing arguments of the scalar function 9. With this we 
can give explicit expressions for Eq. (|273j) . In App. iGl it 
is shown how the conservation of Hermicity by the self- 
energy can be used to minimize the number of terms to 
be calculated. 

For the right hand side of Eq. (|273|) we now give ex- 
plicitly half of terms: 

" \X<r){Xa\ (277) 

[MiiaUP^la^^T+lP^T^ia+^P^la^)) 
+Af((aUP^\a^),(S,\Pf 2 \S a ),(at (T \Pt''\aZ a ))] 
+ |*r)(Xff| (278) 
N ((4 CT |Pf + >+ CT ), ( Xs |P*| X .), (a± g \P™-*\aZ*)) 
-\x*){Zl\ (279) 
[AA((a|jP 1 Q ^|a; (T ),(T + |P 1 T 2 |T + ),(a+ ff |P 2 ^| a + ff )) 
-M {{a% a \P?+°\a- +a ), (S„\P? 2 \S a ), (aljP a a -|ai <r )) 
+M((aUP^\a^)Ax,\P^\x,)A^ s \P^n^)^ 

[N((aX a \P^\al a ),(T + \Pl 2 \T + ),(aX a \P 2 + '\a^)) 

+M((aX a \P^\al„)Ax.\P? 2 \x S ):(aUPt°\^ a ))] 

-|T_)(T_| (280) 

[N{(aUPt°\^ a )ASAP?2\S*),{aUPt°\aZ a )) 

(281) 

+U((a.,\P^\aZ,),(x<r\P 1 x 2 \x S ):(aUP 2 -'\aZ a ))] 
+ (282) 
[AA((T + |P 1 T |T + )(a+jP 1 Q -|a;J,(5 (T |P 3 s |^)) 

+M((T + \PnT + ),(aUP^\a^),( X ,\P 2 x \x,))] 

+ \aZ a )(at a \ (283) 

W {{x„\pf\xA (<*up?r kj, osy^fiso) 

+N({ x ^P?\ Xa ),{a±^ 3 \azMT-\P%\T-)) 
+ M((S a \P?\S a ),(a± a \P^\aZ a ),(T_\Pf\T_)) 

+ Af ((S a \Pf\S a ), (a+jP^KJ, ix S \Pi\x,)) 
-M((x«\P?\x*),(aUPil + *\*^),(x<r\P?\Xt))] 

(284) 

Here we use the notation 

P* = P k (- Vl ^), Pti = P k (-ViVi - Wj) (285) 

We fixed the particle-hole index rj in the multiindices 1 
and 2 to be + (and therefore — for the 1, 2): 



+, ^0') = 1,2 
-, = 1,2 



(286) 



All other indices in Eq. (|277[l - (|284[) are implicit summed 
over. The other half of the terms of Eq. (|273|) can be 
constructed in the same way by taking into account the 
opposite sign of rji using the recipe of Appendix [G]) . For 
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FIG. 6: Zero-temperature non-linear conductance dl/dV 
vs bias V and gate voltage V g = — e for strong interaction 
U = 30F and (a) zero magnetic field B = and (b) finite 
field B = 9F. Both figures are calculated in 2-loop RG ( 
Eq. (|228l l, (12731) ) and are converged with respect to the num- 
ber of non-equibrium Matsubara axes (cf. Sect lHI B lfl using 
a sufficiently large bandwidth D = 10 T. The linear conduc- 
tance in (b) shows good agreement with the Friedel sum rule, 

cf. Fig. m 
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FIG. 7: Linear conductance df/dVlv^o as function of the 
gate voltage V g , obtained numerically for bias V = O.OOir 
for which the response was checked to be linear for all V g . 
Our 2-loop RT-RG results (black curves) are compared with 
Friedel sum rule conductance ^ CT g a (red curves) which are 
obtained from Eq. (|287[) using the occupations (n CT ) calcu- 
lated within our RG. Across panels (a)-(d) the interaction 
increases, U/T = 0.0, 2.5, 5.0, 10.0. Within each panel the 
magnetic field is increased, B/Y = 0.0,1.0,5.0,10.0. Panel 
(a) for U = 0.0 serves as a reference, numerically confirming 
the analytic result Eq. (|253[) that already in 1-loop we include 
the exact non-interacting solution. 



non-interacting limit, {7 = 0, this implies that all 2- 
loop (and higher loops) corrections relevant to the cur- 
rent vanish exactly cf. Eq. (|245[) . We emphasize there 
are additional terms not listed in Eq. (|277|) - (|284j) which 
are irrelevant for the current. These describe the 2-loop 
renormalization of the £ coefficient: these involve fac- 
tors (a~ \P a v" |a+ ) and therefore do not vanish, even 
for C/ = 0. However corrections to this coefficient beyond 
2-loops order do contain the factors l|246p and vanish for 
U = 0. See for more details Rcf. K5(J. 



the calculation of the current only the leftmost reservoir 
index n should not be summed over. 

We presented Eq. (|2"77|) - (|2"54"|) in order to show a 
number of important general properties. First, in the 
fermionic sector (|28ip - (|284[) , the structure is the same as 
for the 1-loop equation, i.e., only terms |a~ cr )(a^" cr | appear 
as discussed Sec. B 3[) . This property holds in any or- 
der of RG and is manifestation of the general properties 
Eq. (fT33l)- (IT34)) . 

Secondly, Eq. (|277 ]) -([2"8"4" P only lists the terms of the 
2-loop RG equations which are relevant for the current. 
The right hand side of all these equations are propor- 
tional to the matrix elements (a+ cr |P Qr,CT |cC CT ) as was an- 
ticipated in Sec. (jlll C 3[) . This is also a general prop- 
erty which holds in any loop order of the RG. In the 



IV. RESULTS 



In this section we perform a detailed numerical investi- 
gation of the zero-temperature 2-loop RT-RG equations 
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where the right hand 

sides are given by Eq. (|228|) and (|273l) , obtaining the cur- 
rent as indicated below Eq. (|21ip . We focus on the de- 
pendence on the interaction U and the magnetic field B 
as function of both the bias and gate voltage. To clearly 
structure the discussion we first summarize the central 
features of the calculated conductance as exemplified in 
Fig. [12a) and (b) for zero and finite magnetic field B, 
respectively, and assess the limits of applicability. 
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A. Overview and limits of applicability 

At zero magnetic field the dominant features in 
Fig. [12a) are the Coulomb blockade diamonds defined 
by lines along which a SET resonance appears. In our 
calculated results, these dl/dV peaks are broadened on 
the scale T due to non-perturbative tunneling processes 
and have a peak height e 2 /h, i.e., the quantum con- 
ductance. Going into either of the Coulomb blockade 
regimes where the charge is quantized to N = 0, 1 and 2, 
the current decays non-exponentially due to higher order 
tunneling (cotunneling and higher order processes). At 
very small bias, however, the conductance shows a pro- 
nounced anomaly, but only in the N = 1 regime where 
the dot has an unpaired spin. We stress from the start 
that this should not be naively identified with the Kondo 
anomaly of the Anderson model: the correct description 
of the Kondo peak requires 3-loop RT-RG corrections 
which are beyond the scope of this wor k 27 ' 53 . To clarify 
in which regimes of voltages our results apply we first 
discuss the linear conductance through the spin-channel 
a: g a — (dI a /dV) v _ Q , in particular the consistency with 
the Friedel sum rule 



fjo 



— sin 2 (7r(n (T )) 



(287) 



In Fig. [7] it is clearly seen that at zero field B = the 
conductance increasingly violates the Friedel sum rule 
between the two SET peaks for larger U, the violation 
becoming maximal at the particle hole symmetry point: 
our result reaches 4e 2 //i instead of 2e 2 /h. At best, in 
this regime our 2-loop approach can be a starting point 
for further 3-loop corrections containing the log-divergent 
Kondo corrections: However, it should be noted that the 
violation is finite, even at zero T: The key observation in 
Fig.[7]is that beyond a magnetic field B ~ T only a renor- 
malized elastic cotunneling background remains and our 
results rapidly becomes consistent with the Friedel sum 
rule. Clearly, at voltages above B ~ Y the 3-loop Kondo 
renormalization is expected to be negligible compared to 
the 1- and 2-loop corrections that we accounted for here. 
In all further analyses, the low bias regime V < Y will 
thus be ignored for magnetic fields B <Y. We do, how- 
ever, show our results in this bias regime for two reasons: 
(i) knowing the behavior of the 2-loop scheme is of in- 
terest as it presents a starting point for future 3-loop 
calculations and (ii) the behavior of the 2-loop approach 
can be compared with that for other methods in this 
regime. We note, e.g., that for U = 2.5Y the violation of 
the Friedel sum rule is still rather modest, even at zero 
field. Our 2-loop RG thus accounts non-pcrturbatively 
for the strong tunneling effects at zero temperature, cov- 
ering the complete finite bias stability diagram, where 
previous perturbative generalized master / kinetic equa- 
tion approache a 29 ! 30 break down. 

Based on the above, we expect that for a finite mag- 
netic field B > r the 2-loop RT-RG calculations reliably 
address transport features, illustrated in Fig. HJb), at 
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FIG. 8: Peak positions of dl/dV: shown is dI 2 /dV 2 in a 
zoom of Fig. El making the zeros of dl /dV stand out as curves 
separating red (positive) and blue (negative) regions. The yel- 
low lines are guides to the eye obtained by accurate extrap- 
olation of the linear parts of the resonance lines (including 
data points outside the figure). These emphasize the change 
of the slope of the linear parts of the resonance positions, in 
addition to the non-linear renormalizations close to the kinks. 
In (a) kinks occur at V ~ and V w ±17, whereas in (b) they 
occur at V « ±B and V = ±U + B. Note that in (b) there 
is no discernible kink at V — ±((7 — B): at this energy there 
is a SET resonance "hitting" the Coulomb blockade diamond 
edge but there is no onset of inelastic cotunneling, in contrast 
to V = ±.B where there is such an excitation. This signals 
the importance of inelastic cotunneling for the appearance of 
such kinks. 



all applied voltages. Clearly, the SET resonance peaks 
have been split due the Zeeman effect. The zero-bias 
anomaly splits into two inelastic cotunneling resonances 
at finite bias V ~ B (Zeeman excitations). A much 
smaller zero-bias anomaly remains which should be ig- 
nored, as mentioned above. The above mentioned fea- 
tures are of course known from previous studies and have 
been observed in many experiments. Our new approach, 
however, includes renormalization effects of these basic 
transport signatures, which are non-perturbative in Y far 
from equilibrium. The following more detailed analysis, 
bearing the above restrictions in mind, indeed reveals 
several low-tcmpcraturc renormalization effects that may 
be of experimental relevance. 
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FIG. 9: (a) Non-linear conductance dl/dV as in Fig.[6jb) but 
for reduced interaction U = 10F and magnetic field B — 3T. 
(b) Derivative of (a), d 2 I/dV 2 , highlighting the renormaliza- 
tion of SET peak position (see caption Fig.[8jb)), in particular 
at the onset of ICT. 



B. Single electron resonance: level renormalization 
and broadening 

1 . Kinks 

Careful inspection reveals that the SET resonance lines 
in fact change their slope when crossing V — and 
V = U. This can already be seen for small T in Fig. [(Ha) 
where we plot the V-derivative of Fig. to follow the 
peak positions, adding a linear extrapolation. The SET 
resonance lines of the inner diamond change in such a 
way that the charge gap hardly renormalizes for U ^> T, 
although diamond distortions are visible by kinks in the 
linear extrapolations. The charge gap can both be deter- 
mined from the height of the diamond (non-linear re- 
sponse) or from the width of the diamond (linear re- 
sponse) and no significant deviation is found in this limit. 

A related effect arises in a magnetic field: the SET 
slope below (above) the ICT threshold is smaller (larger) 



than the slope of the bare resonance line. As a result the 
SET lines now show a kink at finite V = ±B. We note 
that no such kink is seen at V = U — B: this indicates 
that indeed the ICT is responsible for this effect since 
at V = U — B, in contrast to V = B there is no ICT 
excitation. In (a) at V = there is some non-linearity 
around V = T (small V should be ignored, see above) 
which persists in (b) around V — B in a magnetic field. 

Upon increasing T relative to U these effects are en- 
hanced as shown in Fig. IHa). The edges of the N = 1 
Coulomb-blockade regime tend to bend inwards, towards 
the diamond center. Notably, above the onset of COT 
the slope is slightly larger than that of the bare res- 
onance line. We furthermore observe that this also 
leads to the suppression of the excitation | \) — > |2) at 
— = ^ + U in Fig. [9] (see arrow), which is still 
clearly visible in Fig- [HI see arrows in Fig. [SJa) and (b). 

2. Analysis 

It seems not possible to analytically extract a sim- 
ple physical picture explaining the above non-linearities. 
The following analysis aims to indicate why this is the 
case: we trace back at which stage of the 2-loop RG 
scheme the various effects are generated taking the pa- 
rameter set of Fig. H^b) as a starting point. In Fig. [TU] 
we show the conductance calculated both in 1- and 2- 
loop RG and both with and without converging the cal- 
culations with respect to the non-equilibrium Matsubara 
axes. 

Clearly, the different slopes and non-linearities already 
arise in the 1-loop RG: it is visible from Fig. HUT b). Their 
strength correlates with that of the signatures of ICT 
appearing in the stability map in the various approxima- 
tions. However, for this it is crucial that the Matsubara 
axes are accounted for: this is seen by directly comparing 
Fig. HUT a). Fig. [TUTb) and is confirmed by Fig. [TT] where 
we explicitly plot the difference of former two figures. 

Overall, the 2-loop corrections are most pronounced 
along the SET-regime boundaries and the ICT threshold 
as comparison of Fig.[TU|b) and Fig.fTUtd') and the plot of 
their difference in Fig.fTSlshows. Also in 2-loop order, the 
kinks are most pronounced in the Matsubara-converged 
result. We conclude that the 2-loop fluctuation effects re- 
sult in a non-trivial energy dependence of the vertices and 
Liouvillian which shows up in anomalous features of the 
measurable stability diagram, even for such a simple An- 
derson model of a quantum dot. We emphasize that these 
features are not related to renormalization processes that 
generate the Kondo effect (3-loop, not included here) and 
have an effect at V ~ B > T. 



3. Experimental implications 

Having traced the origin of the change of the slopes 
and the non-linearities of the SET resonances we discuss 



36 



l-H 



" . u 



FIG. 10: Comparison of dl /dV calculated (same color scale 
as Fig. [9]) in the 1-loop (a-b) and 2-loop (c-d) approxima- 
tion neglecting all Matsubara axes altogether (a,c) and fully 
converging using 15 frequency axes (b,d). The 2nd column 
shows the corresponding dI 2 /dV 2 maps, allowing the dl/dV 
peak positions to be followed (see caption Fig. EJb)) m P ar ~ 
ticular, the inelastic cotunneling excitation at V = B evolves 
from a step in (a)-(c) into a peak in (d). Clearly, the non- 
equilibrium Matsubara frequency dependence is responsible 
for the "kinks" in the SET resonance (they already appear in 
1-loop with quantitative modifications in 2-loop). 



their relevance to experiments. In fact, kinks in SET 
resonances are often observed in various type of quan- 
tum dot systems^—. Our calculations indicate that 
tunnel-induced renormalization is a possible mechanism 
for their occurrence, but other (e.g. electrostatic mech- 
anism o i ) should not be ruled out in an experimental 
situation. However, for strong coupling it is physically 
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FIG. 11: Effect of Matsubara axes convergence: plotted is 
the converged result minus the result neglecting all Matsubara 
axes for (a) 1-loop RG (Fig. Ugh) -(a)) and for (b) 2-loop RG 
Fig. UOf d) — (c)). Adjacent red and dark blue regions indicate 
that the correction is an S-shaped curve which, when added 
to a peaked curve, results in a shift the peak position. Clearly, 
the Matsubara frequency dependence has an impact on the 
positions of all resonances and should be accounted for fully. 
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FIG. 12: Effect of 2-loop corrections: shown is the 2- 
loop result minus the 1-loop result in (a) without Mat- 
subara frequency dependence (Fig. UOf c) — (Fig. UOf a)") and 
in (b) for converged Matsubara frequency dependence 
(Fig JlOf d - ) — (Fig. llOf b)l Note the positive corrections to mag- 
nitude of the inelastic cotunneling in (b). 



not unexpected that when ICT sets on the level renor- 
malization significantly changes, resulting in such a kink. 

A direct test of this assumption would be to track the 
Coulomb diamond as function of the coupling strength 
r. In Fig. fT3Ta) and (b) we show predictions for the 
evolution of the SET resonance point for two fixed gate 
voltages, one below and one above the ICT threshold, re- 
spectively. The main observation from such a plot is that 
the peaks evolve along curves that are not simply offset 
by a constant bias. This indicates that the renormal- 
ization of SET resonance becomes increasingly nonlinear 
and a kink must develop. (Note that experimentally T 
may change non-lincarly with control voltages but this 
does not spoil the argument.) It is important to prop- 
erly choose the point above the ICT threshold: depending 
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FIG. 13: Distortion of the stability diagram with increasing 
tunnel coupling Y. (a) Zoom in of Fig. [9] of dI 2 /dV 2 (same 
scales and units) high-lighting the gate voltage dependence of 
the SET and ICT bias thresholds generated by tunnel renor- 
malization. Here U = lO.Or*, B = 3.0F* and T* is the refer- 
ence value of r. In panel (b) and (c) we show the evolution 
of the zeros of dl 2 /dV 2 as the tunnel coupling is increased 
from r = 0.3r* to 3.0r* along a fixed gate-voltage cut in (a), 
marked the vertical dashed line V g = — e = 0.65r* for (b) 
and V g = — e = 3.50r* for (c). Also U and B are kept fixed. 
The dashed linear approximation to the renormalized SET 
positions in (c) is copied with a vertical offset to (b) , showing 
that the renormalization is non-uniform in the gate-voltage. 
This signals a distorted stability diagram. The ICT peak has 
a weaker F dependence that is non-monotonic (not visible). 



on the gate voltage position the peak may renormalize 
stronger or weaker than the peak below the threshold. 
We are aware that experimentally such tuning of T with 
gate voltages may lead to other side effects which may 
be hard to distinguish from the effect. Here the differ- 
ent renormalization of the ICT resonance can be of use, 
which is discussed next. 



C. Cotunneling resonance: gap renormalization 
and reduced broadening 

Having discussed the effect of the ICT on the SET 
resonances, we now study the ICT features themselves 
in more detail. Fig. [H] shows how the inelastic cotun- 
neling resonance exhibits a Zeeman splitting with in- 
creasing magnetic field. Despite the zero temperature, 
the width of the inelastic cotunneling feature is finite 
and clearly seen to depend on the magnetic field, and 
thereby on the voltage at which this resonance occurs. 
This is in contrast to high-temperature 2-loop perturba- 
tion theor y 29 ' 30 where this resonance appears as a ther- 
mally broadened feature at the unrenormalized excita- 
tion energy. At larger voltage V = B, the resonance 
width increases, reflecting a decreasing life-time. This 




FIG. 14: (a) Zeeman splitting of the cotunneling resonance 
(change in blue color) and shift of the SET peaks (red peak) . 
U = 30r, V g = (17/2 - U/3)T = f//6 = 5r = -e (b) dl/dV 
traces of (a) for B/T = 2.0,3.0,4.0,6.0,8.0, 10.0 (c) Zoom in 
of the bias-derivative of (a), i.e., dI 2 /dV 2 showing the res- 
onance position. The peak in dl/dV in (a) follows the line 
V — B (marked dotted line starting at (0,0)) The step in 
dl/dV in (a) instead follows a parallel line that is offset by 
a constant magnetic field (in this case « 0.45F, marked by 
the dotted line starting at (10,10)). For lower field, this step 
position starts to deviate from the linear behavior when the 
peak starts to develop and the transition sharpens. 
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cnergy-dcpcndcncc of the width is generated by our 2- 
loop renormalization since the initial Liouvillian of the 
RG flow, Eq. (|136|) . has imaginary parts that are all ~ T. 

Next, as the magnetic field is reduced, but still on the 
order of several times T the differential conductance de- 
velops a pronounced peak on top of the inelastic tunnel- 
ing step: in Fig. Q3Jb) this is signaled by the onset of a 
negative second- derivative of the current (blue) and clear 
from the conductance traces in Fig. [H|c). It is known 
that part of such a peak on top of the well-known in- 
clastic tunneling step^ is due to non-equilibrium occupa- 
tions^!^— that we also fully take into account. The en- 
hanced conductance at the cotunneling resonance is due 
to the 2-loop renormalization, including the frequency de- 
pendence. Only after including both 2-loop corrections 
and converging with respect to the Matsubara frequency 
axes the ICT resonance evolves from a dl/dV step into 
a peak. This was illustrated (for smaller U) in Fig. [TUJ 
(a)-(c) show no ICT peak (since their derivatives (e)-(g) 
have no zero at the ICT threshold V ~ B) in contrast 
to Fig. HUT d) , see also Fig. [13] These effects thus cannot 
be explained by a simpler, yet still accurate effective pic- 
ture, underlining the importance of the general approach 
developed here. 

This enhancement is not related to Kondo-exchangc 
tunneling, which is known to lead to additional logarith- 
mic enhancements^^: their renormalization is not in- 
cluded in our 2-loop calculations and is expected to be of 
limited importance at this large magnetic fields (several 
times r). The qualitative change of the inelastic conduc- 
tance feature from a step to a peak in Fig. RUT b) implies 
some ambiguity in the extraction of the excitation energy 
as either the inflection point of the step (upper right cor- 
ner) or the peak position (lower left corner). 

Finally, close inspection Fig. ITBT a) shows that the ICT 
bias peak position has a weak gate voltage dependence. 
In contrast, in the high-temperature limit this resonance 
lies at the unrenormalized cotunneling excitation V = B. 
Strong tunneling thus leads to an apparent renormal- 
ization in the cotunneling peak position, even for this 
simple model, c.f^. Although our effective Liouvillian 
contains parameters E a and F va that relate to the mag- 
netic field (compare Eq. (|136[) and Eq. (|137|1 ). these pa- 
rameters depend on energy and their values at several 
different energy scales go into our full result since many 
non-equilibrium Matsubara frequency axes must be ac- 
counted for. 



V. OVERVIEW AND OUTLOOK 

We have studied the standard model of an interacting 
quantum dot, modeled as an Anderson impurity, in the 
low-temperature, non-equilibrium limit. We have calcu- 
lated the effective time-evolution kernel L(z) = L + T,(z) 
for the kinetic equation of the reduced density opera- 
tor, where E is the non-trivial self-energy supcroperator. 
In contrast to many previous studies using such a gen- 



eralized / quantum master equation approach, we have 
calculated this time-evolution kernel using the real-time 
renormalization group (RG) . The equations for the effec- 
tive kernel are integrated as function of a cut-off parame- 
ter A: as A is reduced, the renormalized Liouvillian La in 
principle flows to the exact result, La|a = o = L(z). This 
RG calculation of the effective kernel involves a number 
of key elements: 

• Transitions between all Liouville space supervec- 
tors need to be accounted for. This includes those 
elements of L\ and L(z) which in a perturbative 
calculation of the kernel drop out of the calcula- 
tion due to conservation laws (charge, spin, and 
possibly particle- hole symmetry). The reason for 
this is that as one integrates out energy scales, ef- 
fectively higher-order diagrams are included into a 
renormalization of the kernel which describe virtual 
intermediate states which are less restricted by con- 
servation laws. 

• The dependence of the kernel on the real QD fre- 
quency (E) (Laplace variable conjugate to time) is 
important even in the stationary state. During the 
RG flow the renormalized Liouvillian at frequency 
E self-consistently couples to its action on virtual 
intermediate states at frequencies differing from E 
by multiples of the bias voltage fit, — /jlr. We have 
shown that this non-equilibrium effect leads to sig- 
nificant quantitative corrections and may require 
tens of bias-multiples to achieve convergence. 

• The reservoir-frequency (itu) dependence of both 
the kernel and the vertices becomes important 
when going beyond the leading, 1-loop approxima- 
tion, in addition to the QD frequency E. This de- 
pendence is generated on the imaginary frequency 
axis when the reservoirs are integrated out and it 
may cancel 2-loop corrections calculated without 
frequency dependence. 

We have systematically accounted for the leading fre- 
quency correction within a 1 plus 2-loop approximation 
to the exact RT-RG equations and derived an effective 
RG equation for the time-evolution kernel only. This 
includes the relevant vertex renormalization corrections. 
Importantly, the leading frequency correction of the 1- 
loop renormalization of the Liouvillian L\ was found to 
exactly cancel the 2-loop zero frequency term. 

For the non- interacting Anderson problem (U = 0) we 
found that the current is exactly captured already in the 
1-loop approximation without any frequency dependence 
even though the complete solution (i.e., including the 
density operator) is contained only within the 1 plus 2- 
loop approximation. For the strongly interacting case of 
interest the non-linear transport spectrum (dl/dV stabil- 
ity diagrams) was calculated for a wide range of param- 
eters. The different, intrinsic broadening of the single- 
electron tunneling and cotunneling resonances was cap- 
tured, as well as the zero-temperature renormalization of 
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their positions. As emphasized throughout, due to the 
restriction to 1- and 2-loop diagrams, the limited Kondo 
regime cannot be addressed. This regime has been re- 
cently studied in detail using RT-RG approach based on 
the mapping to a Kondo model. This allowed the en- 
tire crossover from weak to strong coupling to be de- 
scribed^. Our study thus provides a starting point for 
a 3-loop analysis of the non-equilibrium Anderson model 
in which the interplay of Kondo spin-fluctuations and 
charge-fluctuations can be described. 

The RT-RG study benefited a lot from an extensive 
reformulation of the underlying real-time perturbation 
theory in terms of vertex superoperators G and G. Al- 
though originally introduced in the RT-RG^I we have 
revealed their full significance as fermionic field super- 
operators that directly generate the Liouville Fock space 
in complete analogy to closed quantum many-body sys- 
tems: 

• Field super operators obey definite anticommuta- 
tion relations, and a simple fluctuation-dissipation 
relation similar to the underlying field operators. 
The Wick theorem in Liouville space at finite tem- 
perature can be obtained algebraically using re- 
lation (|59")) as usual^. In contrast to the path- 
integral formulation^ this was previously not obvi- 
ous^ and complicated by fermionic parity signal. 

• The causal structure of the theory is reflected by 
the vanishing of 2 out of 4 reservoir correlation 
functions (rather that 1 out of 4 as in the Keldysh 
Green function technique). This results in an expo- 
nential reduction of the terms contributing to the 
time-evolution kernel (additional to the reduction 
in the wide-band limit). The remaining contribut- 
ing diagrams are easily identified by their topology. 

We have extended the use of these fermionic field su- 
peroperators to the perturbation theory underlying the 
RT-RG. This resulting causal representation of the per- 
turbation theory has many advantages: 

• Probability conservation of the kernel is manifest 
term-by-term, allowing non-conserving approxima- 
tions to be easily spotted. In addition, other ex- 
act eigenvectors and eigenvalues of the kernel were 
found, that limit the form of the exact effective 
Liouvillian L(z) as expressed in our central re- 
sult Eq. (fTBT)) . 

• Tcrm-by-term diagrammatic evaluation of the 
wide-band limit. Diagrams that vanish in this limit 
can be identified by their topology and the cut- 
off independence of the remaining diagrams can be 
shown from the fundamental fermionic algebra of 
the fields superoperators. This results in a fur- 
ther strong reduction in the number of contributing 
terms as function of the perturbation order. More- 
over, the advantage of working with the complete 
space of the QD states becomes explicit. 



• Finally, the fundamental importance of the infinite 
temperature limit becomes explicit. It defines the 
Liouville Fock-space vacuum and the perturbation 
theory can be explicitly decomposed into infinite 
temperature renormalization effects and the non- 
trivial finite-temperature contributions. This pro- 
vides a natural starting point for the RT-RG, which 
readily suggests itself. 

Besides their application to the RT-RG, we have already 
found useful application of some of these points in pertur- 
bative studies as discussed here in connection to Ref. HH 
and elsewhere^ for adiabatic driving^. 
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Appendix A: Wick theorem for causal field 
superoperators of the reservoirs J q 

Here we show how the algebraic proof of the Wick theo- 
rem^ for standard field operators directly applies to the 
field superoperators J in the causal representation. In 
this proof one considers the average of n reservoir field 
superoperators Jf' defined Eq. (|52"1) . 

X = Tr(Jf...J«V i ), (Al) 

-ft 

(which is nonzero only for even n) and commutes J^" to 
the left side, using that the fields obey the usual anti- 
commutation relations (flyfj) . We first consider the case 
q n = — and make use of zero trace property of the causal 
representation (|62p . This reduces the average to expres- 
sions of the same form 

X =J2(- 1 ) Nk ' n %nX k , n (A2) 
but with averages over n — 2 operators 

X n ,k =Tr (Ji 1 ■■■Jk-l Jk+l ■■■Jn- iPres) (A3) 

weighted with a contraction function 

%,n = -^S kifl 5+,q k (A4) 
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Here Nk,n — n ~ k — 1 is the number of permutations to 
bring J| fc and J% n together. 

For q n = + we proceed in the same way, except that 
when J^" has been permuted to the far left we apply 
the fluctuation-dissipation superoperator identity (|59|) . 
This transforms the expression into that for the q n = — 
case multiplied with the Keldysh distribution function: 
Eq. (|A2I) applies again but with jk,n replaced by 

7fe,„ = tanh(77„w, i /27 1 rn ) 7 fej „ (A5) 

when the multi-index n of J* 1 reads n = rj n a n r n Lij n . 
We thus find the standard Wick recursion relation 

X =J2(-l) N ^ k , n X k , n (A6) 

with 7 fei „ = jk, n S+,q„ +lk, n 5-.q n which by iteration gives 
the Wick theorem Eq. (|60| with contractions (|64 |) - (|63l) . 

Appendix B: The main properties of the causal 
vertex superoperators G q : proof 

a. "Bare" vertices In the (jll B 2|) the main proper- 
ties (|53| and (|57|) the causal superoperators G and G were 
introduced, namely, their anticommutation relations ([55)1 
and their relation by Hcrmitian conjugation in Liouvillc 
space (f57|) . We now give the proof of the latter rela- 
tion. We note that both relations are fundamental as 
they imply a formal correspondence of causal = G,G 
operators with the usual fcrmionic field operators and al- 
low us to develop the "2nd quantization" technique for 
fcrmionic Liouvillc Fock-space. 

Super Hcrmitian conjugation is defined relative to the 
scalar product in Liouvillc space (A\B) = TrA^B where 
A and B are dot operators. To prove Eq. (f57| we first 
notice that the "naive" field superoperators (|40p obey 

(S?f ) f = Sf* (Bl) 

where 1 = r\a and 1 = fja. Superoperators with the 
same Keldysh index p are thus conjugated to each other 
in the usual way: Hermitian conjugation is equivalent to 
inverting the particle-hole index r\. This follows by using 
the cyclic property of the trace: 

(A\9f\B) = TtAU^B = (TrBU^A)* 

= {B\<$+\Ay = {B\{$+) ] \A)\ (B2) 
(A\&f\B) = TiA^Bd'l = (Tre^BU)* = (TrB^Acf)* 

= (B\&--\Ay = (B\(& 1 -) f \Ay, (B3) 

where * denotes complex-conjugation and f denotes 
either usual hermitian-conjugation or super-hermitian- 
conjugation depending on whether it acts on operators or 
superoperators respectively. These superoperators, how- 
ever, have the disadvantage that they satisfy no definite 
fcrmionic or bosonic commutation relations, cf. (|45l) . 



The transformed field superoperators (l47l) that include 
the fcrmion-parity sign do anticommutc (cf. Eq. (|49jl ). 
However, they are not related by super-Hermitian conju- 
gation in the usual way: 

{Q'D^vQl (B4) 

which follows from Eq. (|B1[) and the properties (L n ) = 
L n and [L n , G\] = r)G\ of L n = [n, •]_ 

(gpy = ( P L "^) f = (9?) V" = P L "+X = v g\. 

(B5) 

Finally, the causal field superoperators (f5T|) obtained 
from Q v by a Kcldysh-rotation obey both definite 
commutation relations (|53[) and standard Hcrmitian- 
supcrconjugation relations (|57|) : using p 2 = 1 

(^) t = ^£l> (1+,)/ W ( B6 ) 
p 

= 73E^ 9)/2 ^ = Gf (B7) 

p 

key This clearly demonstrates the fundamental advan- 
tage of the causal representation over the other ones. 

Note that due to the p-dependent superoperators en- 
tering in the definition Eq. (|5ip of our superoperator Clo- 
the latter is either a commutator or anticommutator with 
the fermionic operator d'l depending on the argument on 
which it acts. For example G, hcr A = [d 7 l,A] if A is a 
fcrmionic dot operator (odd in charge, see Eq. (|96|) ) and 
Grj^A — [d 7 l,A} + if A is bosonic (even in charge, see 
Eq. (pTj) 1 ). For the G Vta the opposite holds. The cru- 
cial relation Eq. (|55j) nevertheless holds always for the 
superoperators G: one obtains zero in the first case as a 
trace of a commutator and in the second as a trace of an 
operator which is off-diagonal in charge. 

The analogy to usual field operators extends also to 
the transformation behavior under spin-rotations. For 
usual field operators in Fock space dj., ad s transform as 
irreducible tensor operators (ITOs) of rank 1/2 and in- 
dex a/2. Particle (r/ = +) G q +a and hole (?? = —) field 
superoperators with the appropriate pref actors, aG q _ s , 
are similarly irreducible tensor superoperators (ITSOs) 
of rank 1/2 with index a/2. This applies to all the repre- 
sentations of the field superoperators that we used: the 
same holds for Q\ al &G-a an d 

To prove this we note that for any two superoperators 
A p and B p generated by two operators A and B, in the 
same way as G\ (cf. Eq. (|4"0)) ), 

the commutator of superoperators can be expressed 
in the superoperator corresponding to the commuta- 
tor: [Ap,Bp']_ = p5 p . p {[A,B]^) p . This directly shows 
that the field super operators Sff (|40[) transform in the 
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same way as the field operators d\ under spin-rotations: 
[L s *,Sff]_ = ({S u d 1 ]_) p . The superoperators @7|) and 
(|51|) simply inherit this property since spin and charge 
superoperators commute, [L Si ,L n ] = 0. 

Appendix C: Relation of the causal fermionic 
superoperators G q and Prosen's superoperators 

We note that fermionic superoperators (i.e. possess- 
ing anticommutation Eq. (|53[) and Hcrmitian conjugation 
Eq. (|57|) relations) and a Liouville Fock-space similar to 
ours (Eq. ([ 5T ]) and Eq. (|?T|) - ([551) ) were also constructed 
in Ref. [3l] for the problem of steady state of an open sys- 
tem described by Lindblad equation quadratic in fermion 
operators. It can be checked that superoperators cj from 
that work also satisfy Eq. and thus posses causal 
structure. 

In the work^I the opposite order of construction was 
used: i.e. first a Fock-space was constructed as an or- 
dered product of the Majorana operators (which are lin- 
ear combinations of usual fermion creation and annihila- 
tion operators) and then fermionic superoperators were 
defined on it. The fermionic parity superselection rule 
we use was implicitly taken into account by ordering of 
the Majorana operators and by a specific definition of the 
creation and annihilation superoperators. Careful com- 
parison of the Fock spaces shows that superoperators ct 
from the work^I are related to our G v ^ by the following 
unitary transformation: 

4 = -^(G' +tm + G'_ !m ), j = 2m -1 (CI) 

c ) = A= (<5'-,m ~ G '+,m) > J = 2m 

(relations for the superoperators Cj can be obtained from 
Eq. (|C1[) by super Hcrmitian conjugation). Here m = 
1,2,3,... are the indices numerating fermionic channels 
and G v m are the G\ superoperators renumbered by the 
channel index. For the single level model we have only 
spin channels and: 

C = / ^vA-i m = 1 (r9 ^ 
"> m ~ \ G„, t , m = 2 ^ ' 

For the multilevel model with discrete channels k = 
1,2,3,...: 

G' v>m = G^, m = 2k+(a-l)/2 (C3) 

For the case of infinite and especially continuous number 
of channels it is not directly seen from the definition of 
the work^l how one should construct superoperators Cj . 
Then relation (|C1|) can be used as a definition of the cj , ct 

superoperators, since the superoperators G\,G\ have a 
proper definition also in this limit, see Eq. (|54[) . In con- 
trast to our case the superoperators are neither ITOs 



in spin nor in particle-hole space. This fact complicates 
the group-classification we performed in Sec III C 31 The 
causal structure of the kernel in the representational also 
remains implicit which in general plays a crucial role as 
our analysis shows. Finally, we emphasize the different 
scope of application of field superoperators in our work: 
whereas Prosen's approach was formulated to calculate 
steady states of quadratic effective Louvillians, we here 
extend it to the reservoirs with continuous fields as well 
and to simplify the microscopic derivation of the effective 
Liouvillians for non-quadratic problems. 

Appendix D: Schmutz's fermionic superoperators 

We note that in the Ref. I32H341 alternative field su- 
peroperators a^a^ and a a) a\ were introduced. In our 
notation ([15 ]) - ([17 ]) and ([55 ]) . (|10 J) they read: 

|; T=-\ ( m ) 

They are related to our "intermediate" superoperators 
Eq. (@TJ) as follows: 

(the definition of the work^s contains an additional 
inessential p signs). The use of the additional rj- 

1 — rj 

dependent superoperator p~^~ allows one to compensate 
the inconvenient signs p in the anticommutation rela- 
tions Eq. (|49p and in the Hermitian-conjugation relation 
Eq. (|B4[) . thus restoring the usual fermionic algebra with- 
out the use of the Keldysh-rotation: 

[of ,af 2 ]+ = £i,2<W 2 , (a?) f = a\ (D3) 

The spin and particle-hole group transformations of these 
superoperators coincide with those introduced by us. 
However they do not reveal the important property 
Eq. (f55j) which is crucial in our formulation. Also the 
explicit ^-dependence of the sign-prefactor does not al- 
low one to perform a simple Kcldysh rotation of the 
superoperators to recover this property 

Appendix E: Fermion-parity operator and 
superoperator 

In Sec. (jll C 2p the operator Zr, the right eigen su- 
pervector \Zr) of causal field superoperator G\ (cf. 
Eq. (|86p ). turned out to play an important role. Here 
we discuss its additional properties to further clarify its 
physical meaning. 

b. Fermion parity The causal field superoperators 
(|51l) . written out in terms of field operators, read 

G\' = ^ {di • +?(-!)" • <M-l) n } , (El) 
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where we used (— l) i 



-l) n •(-!)". Note that n 



^2 a n <7 is the occupation operator with n a = d^d a . The 
definition of this representation is based on the fermion- 
parity superselection rule, cf. Sec. (jll B 2p . The explicit 
form (|E1[) makes clear that G^f = G\ has 



z R =\{-i) n 



(E2) 



as right zero eigen supervector: upon substitution the 
two terms in Eq. (|Eip simply cancel since (— l) 2n = 1. 
Noting that (-1)" = Y\e tml " =H(1- 2n a ) we recover 



Z R = |(2n t - l)(2n i - 1) = 2^ - 



(E3) 



the result (|90|) in the main text. Clearly, from Eq. (|E2 



(E4) 

which implies the normalization of the supervector |.Zr). 
The eigenvalue equation ([86]) . Gi\Zr) = [cIi,Zr]+ = 
requires to anticommute with all the fermionic fields: 
Zji is therefore the unique operator (up to normalization 
and a phase) that (anti)commutes with all QD bosonic 
(fermionic) operators in the QD Liouville space, in close 
analogy to Grassmann numbers used in functional inte- 
gral approaches^. Eq. (|E2[) most clearly illustrates the 
physical meaning of the operator Zr: we identify Zr as 
the fermion-parity operator (up to a constant), 



2Zr\ti) = \n), n — even 
1Zr\vl) =—\n), ?i=odd 



(E5) 



Finally, that \Zi) = ^1 is a zero eigenvector of G 1 fol- 
lows from (— l) n <ii(— 1)" = — d\. The eigenvalue equa- 
tion ([55)1 . G\\Zl) — [d\,Zi\- — 0, requires Zl, con- 
sidered as an operator, to commute with the QD field 
operators and therefore with all QD operators, implying 
that indeed Zl is proportional to the unit operator. 

c. Spin and charge rotations By construction the 
two independent operators Zl and Zr transform as 
scalars under both spin- and charge-rotations (cf. Ta- 
ble |T|. They must therefore be related the two scalars 
of these groups (Casimir operators of their Lie- algebras), 
S 2 = EiSf and T 2 = E^ 2 : 



Z, 



hi 



(T 2 + S 2 ), Z R = \{T 2 -S 2 ) (E6) 



This relation follows from the fact that these act as unit 
operators in their respective 2-dimensional subspaces: 
S 2 - 3/4£|<7><<r| and T 2 = 3/4(|0)(0| + |2)(2|). 

d. Fermion-parity superoperator We next consider 
the fermion-parity superoperator, defined naturally by 
the right action of the fermion-parity operator on an op- 
erator (identical results follow for the left action): 



U* = •IZr = •(-!)'" 



(E7) 



This unitary and Hermitian superoperator (U' = 
U,U 2 = 1) transforms the two types of field superop- 
crators into each other, 



UG\U = G\. 



(E8) 



We can thus interchange the role of creation and annihi- 
lation supcropcrators in Liouville Fock-space by a linear 
transformation U. This is similar to the field operators 
d(j and d\ that generate the standard Fock-space: the r\ 
index distinguishing can be inverted by a unitary trans- 
formation 



Z7T (Ty Sy ) J p -iir(Ty-Sy) _ J 



(E9) 



^(L Ty -L Sv ) d 



cf. Eq. (|T06|) and Eq. (fT09l) . Note that Kd va = d n<J as 
well, but this is an anti-unitary transformation. 

This result follows by considering an arbi- 
trary fermionic operator F for which (— \) L F = 
(—l) n F(—l) n = —F. The superoperator U transforms a 
commutator of F with any operator to an anticommuta- 



tor and vise versa: defining L } 



F»±»F: 



UL±Um = F • {-l) 2n ± •(-l) n F(-l) n (E10) 
= F • =p • F = (Ell) 

Commutators with any bosonic operator B remain 
unaffected. Since the superoperator G a is a 
(anti)commutator when acting on a fermionic (bosonic) 
operator and vice versa for G, the superoperator U in- 
terchanges these two: UG\U = G\ and- UG\U = G\. 

e. Multi-orbital Anderson models Finally, we indi- 
cate how the operators \Zi) and \Zr) can be constructed 
for more general multi-orbital Anderson-type models. 
The super vacuum state is 

\Zl) = 2^1, 

where N is the number of orbitals and the prefactor takes 
into account normalization (Zt\Zl) = 1. Eq. (|87|) is then 
simply extended to the maximally occupied state with 
respect to this vacuum 

k=\ a \ r] / 



Using Eq. (|ET2|) and Eq. flU]) Eq. (§0]) generalizes to 



implying Z R = jwl. Here n = J2ka n o- ^ s the total 
dot particle-number operator. All properties of single- 
orbital Zr operator ((anti)commutation relations with 
fermionic (bosonic) operators, transformation properties 
under charge and spin rotation, etc.) hold for the multi- 
orbital case as well. 



Appendix F: Symmetry of the self-energy 

In this Appendix we derive the symmetry for QD 
Eq. (jlOip in contact with the reservoirs from the global 
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symmetries Eq. (f25|) and (f26|) . Quite generally, a quan- 
tity with QD and reservoir contributions, A tot = A + A R , 
is globally conserved when [A tot , H tot ] = 0. The corre- 



sponding Liouville super operators L J 
jtot _ 1 t j icn a j so commu t c: 



[L A ,i R ]=0, 
[i Atot ,L tot ] =0. 



[A* 



] and 

(Fl) 
(F2) 



The commutator of the local part L A = [A, •] with the 
Laplace transformed evolution super operator H(z) = 
Tri(z — i tot ) _1 p R of the reduced density operator (cf. 

fl 

Eq. (|3ip ) must then vanish: 



L A H(z) = L A Tr 



Tr— 

-R Z 



R Z - 

i 



L totr 



TrL 

fl 



(F3) 



using subsequently TrL y 



Tr[# R , 

fl 



.] = 0, Eq. (|F2]1 . 



and finally L" 4 /9 R = which follows from Eq. (f2~Tj) and 



Eq. (jFip . Clearly, this proof applies also for the time- 
evolution superoperator n without any QD-reservoir in- 
teraction, i.e., for L v = 0. Then, using [L^II^)] = 
[Z/ A ,n (z)] = and taking the commutator of L A 
with the Dyson equation that defines the self-energy 
U(z)' = U (z) - Tio(z)£(z)n(z), we find that 



Appendix G: Hermiticy 

/. Hermitian conjugation superoperator The density 
operator is restricted to be invariant under Hermitian 
conjugation in Hilbcrt space. When considering density 
operators as supervectors in Liouville space this Hermi- 
tian conjugation of an operator then corresponds to a 
superoperator that we denote by K, 



K\A) = |At). 



(Gl) 



It is to be distinguished from the Hermitian conjugation 
of a superoperator. K is antilinear and satisfies 



K = 1 



(G2) 



where X is the unit superoperator. Changing the basis 
in Liouville space by the superoperator K we effect an 
antilinear transformation of a superoperator (denoted as 
c-transform in Ref. [26h : S — > KSK. In the time rep- 
resentation the density operator is invariant under this 
transformation: Kp(t) = p(t), implying for the Laplace- 
transformed density matrix Eq. (|32|) Kp(z) = p(—z*). 
Applying K to the kinetic equation (|35|1 we obtain a con- 
jugation relation restricting the kernel E(z): 



This property holds of course also for initial Liouvillian: 



KLK» = [H, 



.tit 



\H, • 



-L». 



The transformation of the fields Sf p follows by applying 
K to Eq. flM]) using (rft.)t = .trf = (K»)d: 

K^fK = <S{ (G4) 

giving with Eq. g7]), jG2|) and KL n K = —L n 

KQ\K = p~ L "^ = (-l) Ln Gl (G5) 

The transformation of the causal field superoperators G q , 
giving Eq. (|120[> in the main text, follows from Eq. (|51|) : 



KG\K = q(—l) 



L" 



G\. 



(G6) 



g. Simplifications using conjugation The transfor- 
mation behavior of the kernel E(z) under conjugation 
K of the basis vectors, KT,(z)K = — £(— 2*) restricts 
the structure of the contributions to S(z) in the renor- 
malized perturbation theory Eq. (|79p . The RG-equation 
for the Liouvillian Eq. (pOH)) . (p24|) . (j267|) have a similar 
structure (since we can eliminate the renormalized ver- 
tices, even in the 2-loop RG approximation, we can re- 
strict our considerations to only the bare vertices G as 
in Eq. (|79|) ). To make use of this, we decompose it into 
conjugate pairs. To illustrate the idea, consider first the 
1-loop approximation to S(z) in perturbation theory: 



S(z)=^S li (z) 



(G7) 



where we now write the sum over 1 explicitly and 
E n {z) = 7(xi)GiniGi = -Kt- n {-z*)K denotes a 
term in which the multiindex 1 has a fixed value and 
Hi = (z — L — xi) . Using this and that 1 is dummy 
summation variable, we can restrict the summation one 
fixed r\ configuration, e.g., r\\ = +, while manifestly pre- 
serving the structure KYi{z)K = — E( — z*): 



(G8) 



The calculation of the supcrmatrix elements is now sim- 
plified: using the notation of Sec. (|III[) and the anti- 
linearity of K, (A\KSK\B) = (K A\S\KB)* we obtain: 



(« 3 |E(z)|«o) = 



(G9) 



s m+ ((k3|Gi|k2)(«i|Gi|ko)(«i|IIi(z)|ko) 



((A- K3 |G 1 | K2 )*( Kl |G I |A' Ko )*( Kl |n 1 (-z*)| K o)*)) 

In our Liouville-Fock basis the matrices representing 
G q are of course real since these are the field opera- 
tors. When the basis supervectors IK3), |ko) correspond 
to diagonal operators (Zi, \a, So, To), the second term in 
Eq. (|G9[) simply relates to the first one since these super- 
vectors are mapped onto themselves by K, cf. Eq. (|115p . 
Eq. (JG9J) simplifies to 



(«a|S(z)|«o) = J2 S ^+( k ^i\^)(ki\Gi\ko) (G10) 



KE(z)K = -E(- 



(G3) 



((«i|ni(«)|/so)-(Ki|ni(-«*)|Ko)*) 
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which is explicitly imaginary at z = iO as it should 
be. Supervectors corresponding to non-diagonal oper- 
ators (S a ,T v ,a") come in pairs related by inversion of 
both indices rj and a. The superoperator K maps these 
pairs onto each other (with a sign change for v = —ar], cf. 
Eq. (|117j) . Noting also that G only has non-zero matrix 
elements for one superkets of each pair, , cf. Eq. (|119j) - 
(|121[) . we see that in (|G8|) cither only the first or second 
term contributes (or neither) when K3 and/or kq is a 
non-diagonal operator. In this way we have effectively 
eliminated the need to evaluate the terms for the rji = — 
using the conjugation relations. This consideration is 
generalized to 2-loop expressions by adding to Eq. (|G8|) 



J2 S 



ril + u V2 + 



(Gil) 



12 



where Eii22~(z) collects all 2-loop terms with a fixed mul- 
tiindices 1 and 2. The same analysis applies to terms in 
the RG equations Eq. ((203) , (EMD , (EUl ■ 

h. Conjugation properties under RG flow In con- 
trast to the charge and spin transformation properties, 
the conjugation property (|G6[) - associated with the fun- 
damental Hermicity of the density operator - is exactly 
preserved under the RG flow. We first note that the ver- 
tices G = G + only determine the initial value of the Li- 
ouvillian L and therefore are not affected by the ensuing 
RG flow: they therefore simply obey Eq. (|G6|) . However, 
the vertices G = G~ flow together with the Liouvillian L 
and thereby acquire a dependence on the dot frequency 
z. These vertices obey the following generalization of 
Eq. (|G6|) to non-zero z: 



KG\{z)K = q {-l) Ln G\{- 



(G12) 



To prove Eq. (|G12|) wc use that the bare vertices, provid- 
ing the initial values of the RG flow, possess the property 
(|G12p . It remains to show that for each scale A the in- 
finitesimal correction to the G generated by the RG flow 
also possesses this property. We therefore apply K • K to 
the RG Eq. (gDH) for vertex G a , insert Eq. i[G2]) . using 
Eq. (|G3[) and assume that the property Eq. (|G12[) holds: 



dA 



k^k= ( _i)c*+i)L"+*<* + i)/2 ( 




^1 * , ¥7 z^G 2 ---G & ---G k _ l 3^"T7 (G13) 




z* +L(-z*) 



(G14) 



n 7(^i) r< r< n n 

G l 1 f? ^-0 2 ...0 Q ...Gfc_l =-; — Crfc 



-z\ - L{-zt) 



(G15) 



At the first equality a sign factor arises when commuting 
of all (— l) L factors to the left, using the fermion-parity 
property (— \) L G = G( — 1) L +1 which is preserved un- 
der the RG as well (since both sides of Eq. f|209[) have 
an odd number of G's). At the second equality we in- 
verted all dummy multiindices i — > i and inverted the 
integration variable u>i ~ ¥ giving a sign (— \) k / 2 due 
to the k/2 antisymmetric contraction functions (counting 
7's and dj/dA), where k is the even number of vertices 
other than G a . Since k/2 and k(k + l)/2 have opposite 
parity for even k the result follows. 
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